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ABSTRACT 

The ISM is governed by supersonic turbulence on a range of scales. We use this simple fact 
to develop a rigorous excursion-set model for the formation, structure, and time evolution of 
dense gas structures (e.g. GMCs, massive clumps and cores). Supersonic turbulence drives the 
density distribution in non-self gravitating regions to a lognormal with dispersion increasing 
with Mach number We generalize this to include scales > h (the disk scale height), and use it 
to construct the statistical properties of the density field smoothed on a scale R. We then com- 
pare conditions for self-gravitating collapse including thermal, turbulent, and rotational (disk 
shear) support (reducing to the Jeans/Toomre criterion on small/large scales). We show that 
this becomes a well-defined barrier crossing problem. As such, an exact "bound object mass 
function" can be derived, from scales of the sonic length to well above the disk Jeans mass. 
This agrees remarkably well with observed GMC mass functions in the MW and other galax- 
ies, with the only inputs being the total mass and size of the galaxies (to normalize the model). 
This explains the cutoff of the mass function and its power-law slope (close to, but slightly 
shallower than, —2). The model also predicts the linewidth-size and size-mass relations of 
clouds and the dependence of residuals from these relations on mean surface density/pressure, 
in excellent agreement with observations. We use this to predict the spatial correlation func- 
tion/clustering of clouds and, by extension, star clusters; these also agree well with obser- 
vations. We predict the size/mass function of "bubbles" or "holes" in the ISM, and show this 
can account for the observed HI hole distribution without requiring any local feedback/heating 
sources. We generalize the model to construct time-dependent "merger/fragmentation trees" 
which can be used to follow cloud evolution and construct semi-analytic models for the ISM, 
GMCs, and star-forming populations. We provide explicit recipes to construct these trees. We 
use a simple example to show that, if clouds are not destroyed in '-^ 1—5 crossing times, 
then all the ISM mass would be trapped in collapsing objects even if the large-scale turbulent 
cascade were maintained. 

Key words: galaxies; formation — star formation; general — galaxies; evolution — galaxies; 
active — cosmology; theory 



1 INTRODUCTION 

The origins and nature of structure in the interstellar medium (ISM) 
and giant molecular clouds (GMCs) represents one of the most im- 
portant unresolved topics in both the study of star formation and 
galaxy formation. In recent years, there have been several major 
advances in our understanding of the relevant processes. It is clear 
that a large fraction of the mass in the ISM is supersonically turbu- 
lent over a wide range of scales, from the sonic length (~ 0. 1 pc) 
through and above the disk scale height (~kpc). A generic conse- 
quence of this super-sonic turbulence ~ so long as it can be main- 
tained - is that the density distribution converges towards a log- 
normal PDF, with a dispersion that scales weakly with mach num- 
ber (e.g.'Vazquez-Semad eni|1994[|Padoan et al.|19"97)|Scalo etaL] 
[1998, Ostriker et al. 19991 

Without continuous energy injection, this turbulence would 
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dissipate in a single crossing time, and the processes that "pump" 
turbulence (generally assumed to be related to feedback from mas- 



sive stars) remain poorly understood (see e.g. 


Mac Low & Klessen| 


'2004[|McKee & Ostriker 2007 ; Hopkins et al. 


2012 and references 



therein). However, provided this turbulence can be maintained, it is 
able to explain the relatively small fraction of mass which collapses 



under the runaway effects of self-gravity and cooling fVazquez- 
Semadeni et al. 2003 ' Li et al.|2004||Lr& Nakam ura 2006; Padom 
& Nordlund 2011). In this picture, star formation occurs within 
dense cores, themselves typically embedded inside giant molecular 
clouds (GMCs), which represent regions where turbulent density 
fluctuations become sufficiently overdense so as to be marginally 
self-gravitating and collapse (Evans 1999, Gap & Solomon|2004[ 
[Bussmann et al.|2 008 1. Some other process such as stellar feedback 
is believed to be responsible for disrupting the clouds, after a few 
crossing times (e.g. [Evans et al.|2009| . The turbulent cascade has 
also been invoked to explain GMC scaling relations, such as the 
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size-mass and linewidth-size relations ( |Larson|1981[|ScovilIe et alT] 
[T987l >. 

However, despite this progress, there remains no rigorous an- 
alytic theory that can simultaneously predict these properties, as 
well as other key observables such as the GMC mass function, and 
the spatial distribution of gas over and under-densities in the ISM. 

The approximately Gaussian distribution of the logarithmic 
density field, though, suggests that considerable progress might be 
made by adapting the excursion set or "extended Press-Schechter" 
formalism. This has proven to be an extremely powerful tool in 
the study of cosmology and galaxy evolution. The seminal work 
by [Press & SchechtCT] ^1974| ) derived the form of the halo mass 
function via a simple (albeit somewhat ad hoc) calculation of the 
mass fraction expected to be above a given threshold for collapse, 
expected in a Gaussian overdensity distribution with the variance 
as a function of scale derived from the density power spectrum. 
|Bond et al.| jl991^ developed a rigorous analytic (and statistical 
Monte Carlo) formulation of this, defining the excursion set for- 
malism for dark matter halos. Famously, this resolved the "cloud 
in cloud" problem, providing a means to calculate whether struc- 
tures were embedded in larger collapsing regions. Since then, ex- 
cursion set models of dark matter have been studied extensively: 
they have been generalized and used to predict - in addition to the 
halo mass function - the spatial distribution/correlation function 
of halos ( |Mo & Wh"itel|1996| >, the distribution of voids (Sheth^ 
|van de W eygaert'2004 1, the evolution and structure of HII regions 
in re-ionization (Haiman et al.^200d{ [Furlanetto et al.|2004| l, and 
many higher-order properties used as cosmological probes. By in- 
corporating the time-dependence of the field, they have been used 
to study the growth and merger histories of halos and to construct 
Monte Carlo "merger trees" ( [Bower|p^9T[ [Lacey & Cole||1993[ l. 
These trees formed the basis for the extensive field of semi-analytic 
models for galaxy formation, in which analytic physical prescrip- 
tions for galaxy evolution are "painted onto" the background halo 
evolution (e.g. [Somerville & Kolatt|[T999l [Cole et al. 2000 ). It is 
not an exaggeration to say that it has proven to be one of the most 
powerful theoretical tools in the study of large scale structure and 
galaxy formation. 

There have been other, growing suggestions of similarities 
between the mathematical structure of the ISM and that invoked 
in excursion set theory. The mass function of GMCs, for exam- 
ple, has a faint-end slope quite similar to that of galaxy halos 
(both close to dii/dM oc M~^), suggestive of hierarchical col- 
lapse. [^zquezjSemadeni]|T994|l attempted to rigorously examine 
whether the structure of the ISM should be "hierarchical," although 
they strictly define this as the probability of many independent fluc- 
tuations dominating the "peaks" in the density distribution (which 
does not technically need to be satisfied in excursion set theory). 
This is related to (but not equivalent to) the large body of work 
on the quasi-fractal structure of the ISM (see e.g. 'Elmegreen'2002' 
and references therein). On smaller scales, Krumholz & McKee 
(j2005j suggested that the fraction of a lognormal PDF above a "col- 
lapse threshold" at the sonic length could explain the fractional 
mass forming stars per free-fall time, inside of GMCs. |Padoan| 
|et al.| \\991) \ [Padoan & Nordlund| suggested that the dis- 

tribution of lognormal density fluctuations above a threshold over- 
density could explain the shape of the stellar initial mass func- 
tion (IMF). [Scalo et al.|jI998^ explicitly discuss the analogy be- 
tween this and cosmological density fluctuations, an d|Heiuiebelle| 
|& Chabrier| ( [2008] > expanded upon the [Padoan et al.| ( |1997 l| argu- 
ment using a derivation almost exactly analogous to the original 



[Press & Schechter] ( |I974[ > derivation, and showed that it agreed well 
with the standard IMF. 

But despite these suggestions, and the enormous successes 
of the excursion set model in cosmological applications, there has 
been no attempt to translate the excursion set formalism to the prob- 
lem of the ISM and GMC evolution. At first glance, it is obvious 
why. The cosmological excursion set theory is applied to small fluc- 
tuations of the linear density field, in the linear regime, to dark mat- 
ter (coUisionless) systems, with Gaussian, nearly scale-free fluctu- 
ations seeded by inflation, and to Lagrangian "halos" which (mod- 
ulo mergers) are conserved in time. The Gaussian distribution of 
ISM densities represent large fluctuations in the logarithmic density 
field, which are a product of a fully non-linear, turbulent, gaseous 
(collisional) medium, and evolve both rapidly and stochastically in 
time. 

However, in this paper, we will show that although the physics 
involved are very different, none of these differences fundamentally 
invalidates the underlying mathematical formalism of the excursion 
set theory. 

Here, we develop a rigorous excursion-set model for the for- 
mation, structure, and time evolution of structures in the ISM and 
within GMCs. We show that this is possible, and that it allows us to 
develop statistical predictions of ISM properties in a manner anal- 
ogous to the predictions for the halo mass function. In §[2] we de- 
scribe the model. First (§ [2.1^ , we derive the conditions for self- 
gravitating collapse in a turbulent medium (the "collapse thresh- 
old"), in a manner generalized to both small (sonic length) and 
large (above the disk scale height) scales. Next (§ [2.2^ , we dis- 
cuss the density field, and, assuming it has a lognormal character, 
construct the statistical properties of the field smoothed on a phys- 
ical scale R, which allows us to define the excursion set "barrier 
crossing" problem. In § |3] we use this to derive an exact "self- 
gravitating object" mass function, over the entire range of masses 
(from the sonic length to disk mass), and show that it agrees re- 
markably well with observed GMC mass functions and depends 
only very weakly on the exact turbulent properties of the medium 
(including deviations from a lognormal PDF). In §|4] we show that 
the model also predicts the linewidth-size and size-mass relations 
of GMCs, and their dependence on external galaxy properties. We 
also examine how this depends on the exact properties of the tur- 
bulent cascade. In § |5] we extend the model to predict the spa- 
tial correlation function and clustering properties of clouds (and, 
by extension, young star clusters), and compare this to observa- 
tions. In §|6] we predict the size and mass distributions of under- 
dense "bubbles" or "holes" in the ISM which result simply from the 
same normal turbulent motions. We show that this can explain most 
or all of the distribution of HI "holes" observed in nearby galax- 
ies, without explicitly requiring any feedback mechanism to power 
the hole expansion. In § |7] we generalize the model to construct 
time-dependent "GMC merger/fragmentation trees" which follow 
the time evolution, growth histories, fragmentation, and mergers of 
clouds. In § [7.2[ we provide simple recipes to construct these trees, 
and discuss how they can be used to build semi-analytic models 
for GMC and ISM evolution and star formation, in direct analogy 
to semi-analytic models for galaxy formation. We use a very sim- 
ple example of this to predict the rate at which the gas in the ISM 
collapses (absent feedback) into bound structures, show that this 
agrees well with the results of fully non-linear turbulent box sim- 
ulations, and argue that feedback must destroy clouds on a short 
timescale (a few crossing times) to prevent runaway gas consump- 
tion. Finally, in §[8] we summarize our results and conclusions and 
discuss a number of possibilities for future work, both to improve 
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the accuracy of these models and to enable predictions for addi- 
tional properties of the ISM. 



2 THE MODEL 

The fundamental assumption of our model is that non-rotational 
velocities are dominated by super-sonic turbulence (down to some 
sonic length), with some power spectrum P{k) or which is 

maintained by any process (presumably stellar feedback) in ap- 
proximate statistical steady-state. As we discuss in § [8] all other 
assumptions we make are convenient approximations to simplify 
our calculations, but it is possible to generalize the model. 

The two key quantities we need to calculate the cloud mass 
function and other properties are the conditions for "collapse" of a 
cloud (i.e. conditions under which self-gravity can overcome tur- 
bulent forcing) and the power spectrum of density fluctuations. Be- 
low, we show how these can be calculated for a turbulent medium 
from the velocity power spectrum; however, in principle they can be 
completely arbitrary (for example, specified ad hoc from numerical 
simulations or observations). So long as they are known, the rest of 
our model proceeds identically. 



2.1 Collapse in a T\irbulent Medium 

First, for simplicity, consider gas in a galaxy whose average prop- 
erties are evaluated on a scale R where the velocity dispersion is 
highly supersonic (R ^ ^sonic, where £sonic is the sonic length), but 
where shear from the disk rotation and large-scale density gradients 
can be neglected (R <^ h, where h is the disk scale-height). The tur- 
bulent dispersion on these scales is (v, (A)^)~yt£(yt). If the turbu- 
lence has a power-law cascade over this interval then E{k) oc k^''. 
If the region has some mean density p (on the same scale R), then 
the potential from self-gravity is \U\ ^ jSGM/R « l3Gp{AiT /3)R- 
while the kinetic energy in turbulence is (1/2) (v,(/?)") ; the re- 
gion will be gravitationally bound and "trapped" when pc ~ 
{3 /Sir 13) {v,{Rf) /GR^, where 13 ~ 1 depends on the shape (inter- 
nal structure) of the density perturbation. Formally, we also need to 
check whether the momentum "input" rate from the turbulent cas- 
cade (equal to the dissipation rate in steady-state) is less than the 
gravitational force, and whether the energy input rate is less than 
the rate at which a gravitationally collapsing object will dissipate. 
However, because for super-sonic turbulence, the timescale for en- 
ergy or momentum dissipation on a scale R just scales with the 
crossing time fcoss = R/vt(R), we obtain the identical dimensional 
scaling for all of these criteria|^ 

These are simply a restatement of the Jeans criterion, for 

' There are different conventions in the turbulence and excursion-set liter- 
ature for the normalization and i:-dependence in the definition of P{k). To 
simplify matters, we will refer to the velocity power spectrum by means of 
E{k), which with the assumption of isotropic turbulence gives the differen- 
tial energy per mode as d£ = E{k) dk. 

^ The energy injection rate in the turbulence is = (1/2) i'^(S)/ (77 /cross) = 
{l/2)vf /{riR/vt) where ~ 1 is constant. A virialized object where 
cooling is rapid (i.e. pressure forces can be neglected), where the 
virial motions are turbulent, will then just lose energy at a rate = 
(1/2) \U\/{riR/^l3GM/R) - equating these gives an identical dimen- 
sional requirement on p to the binding criterion, but with a slightly differ- 
ent coefficient. Equating the turbulent momentum input rate d(Miv)/df = 
Mv,/(r; (cross) to the gravitational force Fg,nv ~ GM/R^ again gives the 
identical result. We should take the most stringent normalization from these 
as the relevant criterion, but this is entirely degenerate with the value of 



wavenumber k^ l/R, but with the sound speed Cj replaced by the 
turbulent velocity dispersion v,. For an individual fc-mode (sinu- 
soidal density perturbation), the criteria becomes 



4nG 



oc k^~'' oc R 



p-3 



(1) 



where the latter equalities assume a power-law spectrum (Vazquez^ 
[Semadeni & Gazol 19951. If the system is marginally stable 
with density po on scale Ro, then this simply becomes p{R) > 
po {R/Ro)''~^- If we are in the super-sonic regime, then we expect 
something like Burgers turbulence ( |Burgers|1973[ l, with p ~ 2; but 
we will discuss this further below. 

Now generalize this to a more broad range of radii. On small 
scales, we need to include the effects of thermal pressure: this 

amounts to a straightforward modification of the Jeans criterion 

,.2 



with vf ci + vi jChandrasekhar|l95l|[Bonazzola et al.|l 987l)p] 
On large scales, we need to include the effects of rotation stabiliz- 
ing perturbations. If we focused only on very large {R 2> h) scales, 
where we can neglect the disk thickness, then we simply re-derive 
the ( iToomre|1977^ dispersion relation and collapse conditions, with 



the gas "dispersion" (jg = v? + cJ. More generally, Begelman & 



IShlosman] (12009 ) note that the dispersion relation for growth of 
density perturbations in a turbulent disk (with finite thickness h) 
can be written: 



(2) 
(3) 



where S = 2/ip is the disk surface density, p the average density 
on scale k, and h the disk scale-height, v, the turbulent velocity 
dispersion, and k the usual epicyclic frequency. This differs from 
the infinitely thin-disk dispersion relation by the term (1 + |fc| /?)"', 
w hich accounts for the finite scale height for modes wit h scales A < 



Romeo|l992] nNote that 



h ( |Vandervoort[1970[ 'Elmegreen 1987 
this relation nicely interpolates between the Jeans criterion which 
we derived above on small scales (k ^ h~'), and the Toomre (thin- 
disk) dispersion relation on large scales (k <C A~'). 

If the average density is po, and corresponding average surface 
density Eo, then we can define the usual Toomre Q at the scale h 



Q0{h) : 



ttGEo 



ttGEo 



(4) 



where the second equality follows from (Jg{h) = hfl, which is be 
true for any disk in vertical equilibrium, and we define k = k/^ 
(= y/2 for a constant-Vc disk). If we define the convenient dimen- 
sionless form of k,k=\k\ h, we can write the criterion for instability 
(oj- < 0) as 



po 



> 



po 



Qojh) 
2k 



(1+^)[J 

la. 



k + Kk" 



(5) 



/3 ~ 1. For a rigorous derivation of each of these criteria, see |Bonazzola| 

[snq{T987) . 

it is likely that the power spectrum of velocities Vt will change as we go to 
scales below the sonic length; however, since (by definition) i', < c, in this 
regime, such corrections have essentially no effect on our results. Moreover 
the change - expected to be e.g. a transition from p = 2 to p = 5/3, is small 
for our purposes. 

* Eqn.plis an exact solution for a disk with an exponential vertical profile. 
It is also always asymptotically exact at small and large \k\ and tends to 
be within ^ 10% of the exact solution at all \k\ for the range of observed 
vertical profiles (see |Kim et al.|20Q2^ . 
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Note that the assumption of a finite Qq ensures that so long as there 
is any non-gaseous component of the potential, the gas alone is 
not self-gravitating on arbitrarily large scales (this is important be- 
low, to un-ambiguously define the largest self-gravitating scales of 
clouds). Again, on small scales kh^ I, this reduces to the Jeans 
criterion pc — (jg(k)~ / {A-k G) oc R''~^ , and on large scales kh <^ 1 
it becomes pc = Ec/2h = (kh)~^ / (AttG) oc R. 

|Kim et aL| l |2002| note that is straightforward to further gener- 
alize this criterion to include the effects of magnetic fields by taking 
CTj = V? + Cj + Va, where is the Alfven speed. If we follow the 
usual convention in the literature and assume fi = c]/v\ is con- 
stant, then changing the strength of magnetic fields is identically 
equivalent to changing the sound speed/mach number (which we 
explicitly consider below). Even if we allow /3 to have an arbitrary 
power spectrum, the results are quite similar to this renormaliza- 
tion - for any power spectrum where the magnetic energy density 
is peaked on large scales, it is nearly equivalent to renormalizing the 
turbulent velocities; for a power spectrum peaked on small scales, 
equivalent to renormalizing the sound speed. We therefore will not 
explicitly consider magnetic fields in what follows, but emphasize 
that they are straightforward to include if their power spectrum is 
known. 

Formally, the turbulent velocity power spectrum E(k) must 
eventually flatten/turn over on large scales R^h, both by defini- 
tion (since h itself traces the maximal three-dimensional disper- 
sions) and to avoid energy divergences. If it did not, we would 
recover v, > Vc on large scales in gas-rich systems! Constancy 
of energy transfer and energy conservation require that the slope 
become at least as shallow as E{k) cc k~K A good approxima- 
tion to the behavior seen in simulations is obtained by general- 
izing the exact correction for k near the lowest wavenumbers in 
the inertial scale in Kolmogorov turbulence jBowman|1996| >, tak- 
ing E{k) ^ E{k) (1 + |yt/j|"^)''"'''/^, which interpolates between 
these regimes. This may not be exact. Fortunately however, even if 
we ignored this correction entirely, we can see immediately from 
Equation |5] that for any reasonable power spectrum {p < 3), the 
dominant velocity/pressure term on scales > h is the disk shear 
(~ kR), not Vf. We therefore include this turnover, but stress that it 
is not necessary to our derivation and has only weak effects on our 
conclusions. 



2.2 The Density Distribution 

The other required ingredient for our model is an estimate of the 
density PDF/power spectrum. We emphasize that the our method- 
ology is robust to the choice of an arbitrary PDF and/or power 
spectrum in p. We could, for example, simply extract a density 
power spectrum (or fit to it) from simulations or observations. This 
is, however, less predictive - so in this paper, we chose to focus on 
the case of supersonic turbulence in which case it is possible to (at 
least approximately) construct the density PDF knowing only the 
velocity power spectrum information. 

As discussed in § [T| in idealized simulations of supersonic 
turbulence with a well-defined mean density po and mach number 
on a scale ~ l/R, the distribution of densities tends towards a 
lognormal distribution 



where because po is the mean density. 



dp{5\k) = 



1 



po 



exp 



(--) 



:ln 



-In 



d(5 



(6) 
(7) 



In 



-I 



(8) 



This form of the PDF and our results are identical whether we de- 
fine all quantities as volume-weighted or mass-weighted, so long 
as we are consistent throughout: here it is convenient to define all 
properties as volume-weighted (otherwise po is scale-dependent). 

The dispersion in these simulations is a function of the 
rms (one-dimensional) Mach number averaged on the same scale 
M(k). 



<7,R. ( ln[l + ^M(A:)^ 



1/2 



(9) 



which is naturally expected for supersonic turbulence with efficient 
cooling (because the variance in In (p) in "events" - namely strong 
shocks - scales as In (A^^))|^ 

If the turbulence obeys locality - i.e. if the density distribu- 
tion averaged on some small scale R\ depends only on the local 
gas properties on that scale as opposed to e.g. the structure on 
much larger scales R2 ^ Ri ~ then the distribution of densities 
S{x, R) averaged over any spatial scale R with some window func- 
tion IV (x, R) is also a lognormal in S, with variance 



{R)= d\n{k)al{M[k])\W(k,R)\ 



(10) 



where W{k, R) is the Fourier transform of W{x, R). This is easy to 
see if we recursively divide an initially large volume (e.g. the en- 
tire disk) into sub-regions with different mean po and A4 on scale 
R; each of these sub-regions is a "box" that should obey the den- 
sity distributions above, and so on. Because it greatly simplifies the 
algebra, we will generally follow the standard practice in the excur- 
sion set literature and choose W{k, R) to be a Fourier-space top-hat: 
W{k\R„) = lifk< and W{k \ = if yt > i?,;'. This choice 
is arbitrary, but so long as it is treated consistently, our subsequent 
results are essentially identical (we will show, for example, that us- 
ing a Gaussian window function makes a small difference in all 
predicted quantities)]^ 

It should immediately be clear, however, that if we simply ex- 
trapolated A4^ = v^/cj oc the dispersion would be diver- 
gent! Physically, this would imply ever larger fluctuations in logp 
on arbitrarily large scales; but this cannot be true once the scale R 
approaches that of the entire disk. As kh — )■ 0, the fact that the disk 
has finite mass means that aii — >■ 0. The resolution of this apparent 
dilemma is evident in Equation|2] what matters in A4 in the disper- 
sion is the effective "pressure" from c^; on sufficiently large scales 
kh <^ 1, the differential rotation K./k plays an identical physical 



^ The exact coefficient in front of Ai^ in this scaling does depend on e.g. 
the form of turbulent forcing and other details iFederrath et al. 2010, Price] 
|eraL][20TTt . For our purposes, however, this is entirely degenerate with 
the normalization of the velocity/scale height of the disk and enters very 
weakly (sub-logarithmically). It is potentially more important, however, on 
small scales near the sonic length. 

* As has been discussed extensively in the EPS literature, this does intro- 
duce some ambiguity in the definition of "mass" in the mass function, since 
the real-space window volume is not well-defined. In practice, if we adopt a 
fixed definition of volume = (47r/3)R'i',,, the corresponding systematic dif- 
ferences are relatively small (< 10%) between different window function 
crossing distributions (see |Zentner|2007| . 
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role. We can therefore generalize Equation|9]as 

3 vj{k) l\i/2 



In 1 



In 1 



4 c? + k2 k- 
3 M^{k) 



'\kh\ 



1/2 



(11) 



(12) 



where Mli, = o'gQi) /cj. This ensures the correct physical behavior, 
(j{ — >■ as A: — )■ for all plausible turbulent E{k). 

At some level, our assumptions must break down. And al- 
though it is well-established that the density PDF at the resolu- 
tion limit in numerical simulations (in a "box-averaged" sense) 
approaches the behavior of Eqns. |6|9[ it is less clear whether we 
can assume this on a k-by-k basis and so derive Eqns. [TO)ll| The 
lognormal character of the density distribution holding on various 
smoothing scales as we assume is, however, supported in the in- 
vestigation s of|Lemaster & Stone|p009^ and _Passot & Vazquez-| 
|Semadeni|fr998> ; |Scalo et al.| ( |1998^ . And any distribution which 
is lognormal in either real space or k space must be lognormal 
in both. Moreover, the robustness of this assumption is supported 
by the conservation of lognormality in resolution studies, since all 
simulations essentially measure the PDF smoothed over a window 
function corresponding to their resolution limits. To the extent that 
there is some violation of these assumptions in e.g. the higher- 
order- structure functions (although they are largely consistent with 
locality when Mi, is large; see Boldyrev et al.|2002; Pad oan et al.| 
|2004[|Schmidt et al.]2008^ , this is really a question of the degree to 
which the density PDF globally departs from a lognormal, which 
we discuss below. 

What is somewhat less clear is how Eq. |9] generalizes on a 
scale-by-scale basis. Analytically, the same arguments that prove 
the density distribution of isothermal turbulence should converge 
to a lognormal with real-space variance = ln(l + (3/4)(A1^)) 
trivially generalize to a k-by-k basis (Eq.|9] see |Passot & Vazquez-| 
[Seinadeni 1998"||Nordlund"&Padoan 1999 1. If locality also holds, 
Eq. |lO|must follow. This is the origin of the expectation for an- 
alytic models of the density power spectrum. Note that, as de- 
fined, erf is equivalent to the logarithmic density power spec- 
trum, = kE[ap{k). When M is not large, erf in Eq. [9] scales 
oc M(ky oc V? ~ kE{k), so Ei„p{k) oc E{k). This is just the well- 
known expectation that in the weakly compressible regime, the log 
density power spectrum should have the same shape as the veloc- 
ity power spectrum. | Kowal et al.| ( |2007] l and |Schmidt et al.]j2009^ 
show that this is a good approximation for the In (p) field in simu- 
lations of supersonic turbulence. At higher A4 , this should flatten 
logarithmically, and this is seen in numerical simulations in |Kowal| 
|et al.|(2007| ), in excellent agreement with Eq.|9] These behaviors, 
and the approximate normality of ln(p), appear to hold even in 
simulations which include explicitly non-local effects such as mag- 
netic fields, self-gravity (excluding the collapsing regions), radia- 
tion pressure, photoionization, and non-isothermal gas with real- 
istic heating/cooling (see e.g. |Ostriker et al.ll999{ |Klessen|2000| 



[Lemaster & St one 2009)|Hopkins et al.|2012|l 

Even if our analytic derivation is not exact, we can think 
of the resulting <y^{R) and implied log density power spectrum 
{E\np{k) ~ k~^al) as a convenient approximation for the power 
spectrum measured in hydrodynamic simulations and observations. 
At sufficiently large k, where M is small, E]„p{k) cck^^M^ ock^''; 
a steep falloff with k for typical p r; 2; at smaller k (but still smaller 
scales than the disk scale-height) M is large and this flattens to 
Einp{k) oc k^'lnM^ oc fc^' with a small logarithmic correction. 
This is exactly the behavior directly measured in numerical simu- 
lations ( [Kowal et al.|2007|[Schmidt et al.|2009| ). QuaUtatively sim- 



ilar behavior is seen in the linear density spectrum, but it is impor- 
tant to distinguish the two, since it is well known that large fluc- 
tuations at higher A4 will further flatten the linear spectrum (see 
Scalo et al. 1998; Vazquez-Semadeni & Garcia 2001 [[Kim & Ryu] 



2005[|Kritsuk et al.|2007t|Boumaud et al.,2010| l. It is also consis- 



tent with observations of the projected surface density power spec- 
trum in local galaxies and star-forming regions ( [Stanimirovic et af] 
[l999[ [Padoan et al.|2006, Block et al. 2010| l. If we integrate to get 
cr(if) we obtain a —^constant as if 0, with an absolute value of 
a{R)^l .25 - 1 .9 dex for a range p = 5/3 - 2 and M;, = 10 - 50. 
This range is quite similar to the range measured in a{R) on the 
smallest resolved scales in a wide range of simulations that have a 
sufficiently large dynamic range in scales to probe the typical Mach 



numbers in GMCs and disk scale heights (see Vazquez-Semadeni 



Igg^i 'Nordlund & Padoan T 999 ; Ostriker et al. 20011 |Mac Low" 
klessen 2004; Slyz et al. 2005 , Hopkins et al. 2012). It also agrees 
well with measured values of the dispersion in the real ISM ( |Wong[ 
[et al.|2008[[Goodman et al.|2009a[[Federrath et al.|20I0p . 



3 THE MASS FUNCTION 

The question of the mass collapsed on different scales is now a 
well-posed barrier crossing problem. The quantity S{R) - the log- 
arithm of the density smoothed on the scale R - is distributed as a 
Gaussian random field with variance (R) and zero mean, with a 
well-defined barrier 



5r{R) = 5{p„R)^\n 



Po 



-In 



(13) 



which, upon crossing, leads to collapse. The mass of a collapsed 
object is simply the integral of the density over the effective vol- 
ume of a window of effective radius R-,,. in real space. If the 
medium were infinite and homogenous, this would just be M{Ru) = 
(47r/3) Pc{Rw)Rw', however, we need to account for the finite verti- 
cal thickness of the disk. For the same vertical exponential profile 
that gives rise to the dispersion relation in Equation|2] the total mass 
inside if,,, is 



M{p\R„) = 47vp{R„.)h' 



+ 1 



Ru 



exp 



Ru- 
~h 



0->l 

(14) 

where p{Rw) is the midplane density (chosen for consistency with 
the dispersion relation). This formula simply interpolates between 

M={4tt/3)pR^ foiR<handM = -K(2ph)R- = nT.R- for R>h, 
as it should. 

The fraction of the total mass which is in collapsed objects, 
averaged over a given smoothing scale R,,., is then just 



FcoU (Rw) 



M{p\R,,)p{S\R,,.)dS 



(15) 



where p{3) = po exp (5 — CT[if„.]^/2). Naively, we would equate 
this to the mass function of such objects with the relation 
Mot dFcoii /AM — MAN (M) /dM. Indeed, up to a normalization fac- 
tor, that is exactly the original approach of [Press & Schechter[ 
l [1974^ . However, this neglects the "cloud in cloud" problem: 
namely, it does not resolve whether or not a collapsed region on 
a scale if 1 is independent, or is simply a random sub-region of a 
larger object collapsed on a scale ifo > if 1 . For the case of a con- 
stant 5c, accounting for this amounts to a simple re-normalization; 
but there is no simple closed-form analytic solution for the compli- 
cated 5c here, and we will show that accounting for this behavior is 
critical. 
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Figure 1. Predicted & observed (orange) GMC mass functions. The gen- 
erally predicted mass function is dimensionless; we normalize it to the 
observed surface density Sgas, gas density (or scale height) no, and total 
gas mass Mgas- Together with the assumption that Q ~ 1, this completely 
specifies the model. For each case, we show the exact (Monte Carlo) mass 
function (solid black), and the mass function if we ignore the "cloud in 
cloud" problem by counting bound mass on all scales (dotted red); and 
the analytic fit to the mass function in Eg. |23| (dashed blue). Top: Milky 
Way. The observed MFs are taken from Williams & McKee| (l997| (solid 
line) and Rosolowsky ( 2005 1 (orange points in each panel), & model nor- 
malized to (Egiis, no-'Wgas) = ( IOMqp c-^, Icm'^. 3 X IO'Mq). Mid- 
dle: LMC. Observed MF from |Fukuie t al. 1 2008 1 (line), normalized to 
(8M0pc"-, 0.8cm"3^ 3 X 10^ Mq) (see ^Wong et al. 2009^1. Bottom: M33 
Normalized to (5 M0 pc~^ " " ' ■- • — 

[2003) . 



3.1 Exact Solution 

To derive the exact mass function solution we turn to the standard 
Monte Carlo excursion-set approach. Consider the density field 
at some arbitrary location x, smoothed over some window corre- 
sponding to the radius R (and mass M) 5(xji?„,). This is the con- 
volution <5(x|i?„) = J d'x'W(lx' — x|,/?,v)(5(x'); so if we Fourier 
transform, we obtain 5(k|/?„) = lV(k | (5(k). In other words. 



1.5 cm 1 X IO'Mq) (see Engargiola et al 



Standard (p = 2, M„ = 30) 

p = 5/3 

Mh = 10 

Mj = 1000 

No Low-K E(k) Cutoff 

Power-Law PDF (dM/dlnp-p 

Gaussian Window Function 




t 0.1 r 



0.01 r 



-1 

Figure 2. Variation in the predicted GMC mass function with model as- 
sumptions. The MFs are plotted in dimensionless units. We compare the 
standard model (from Fig.[TJ, which assumes a turbulent spectral index p = 
2, and Mach number at scale ^ /) of A^;, = 30. Assuming p = 5/3 instead 
slightly flattens the slope at intermediate masses. Changing Aii, = 10, 1000 
increases/decreases the sonic length, below which the MF flattens, but near 
the MF break, the assumption of Q ^ 1 means that Mi, factors out. Re- 
moving the assumed cutoff in the turbulent power spectrum at scales S> h 
makes the cutoff in the MF shallower at large masses. Using a Gaussian 
window function to smooth the density field (instead of the usual i:-space 
tophat) makes the MF slightly more shallow, because for the same window 
volume (same mass definition), the radii which contribute fluctuations are 
shifted. In all these models, the density PDF is assumed to be lognomial; 
if we instead assume it is a pure power-law distribution (Eg. |32) , but as- 
sume the same variance in Inp, the result is nearly identical. In all cases, 
the variations in the MF are very small - the marginal stability assumption, 
and weak (logarithmic) mnning of density variance with scale mean the MF 
shape is lai'gely independent of even substantial model assumptions. 



the amplitude 5(x|i?„ ) is simply the integral of the contribution 
from all Fourier modes 5{k), weighted by the Fourier-space win- 
dow function. 

In this sense, we can think of the (statistical) evaluation of 
the density field as the results of a "random walk" through Fourier 
space. Bond et al.]( |1991[ > show that this integration becomes partic- 
ularly simple for the case of a Gaussian field with a Fourier-space 
tophat window, in which case the probability of a transition from 
Si to $2 = S[+ AS as we step from a scale ki to k2 is given by 



p(Si+AS)dAS ■- 



1 



exp 



{AST 
2AS 



d{AS) 



AS = Si -Si =a-{R2)-a{R,) 



where we define the variance 



S{R)^a\R) 



(16) 
(17) 

(18) 



i.e. the increment A 5 is a Gaussian random variable with standard 
deviation \/A5. 

If we begin on some sufficiently large initial scale k ^ 
(R — >■ <x), then the overdensity S and density variations ct{R) must 
go to zero. We then have the well-defined initial conditions for the 
walk, (5(i?max ^ oo) — 0, S{Rnr,,x ^ oo) = 0. Starting at some arbi- 
trarily large i?max, and moving to progressively smaller scales with 
increment^in R or S (ARi or A5,) we can then compute the tra- 

^ The walks defined in this way will always converge as Ai? — > 0. In prac- 
tice, the value of AR should be sufficiently small to ensure multiple barrier 
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jectory S{R) or 5{S), 



(19) 



At each scale Rj, we then evaluate whether or not the barrier has 
been crossed, 



5{Ri) > S,{Ri) 



(20) 



If this is satisfied, we then associate that trajectory with a collapse 
on the scale Rj and mass M{pc[Ri] \ Ri) = M{Ri). 

Recall, we are sampling the field 5(x|iJ„,), so the fraction of 
trajectories that cross the barrier in some interval Ai?, or (equiv- 
alently) AM(^,) represents the probability of an Eulerian volume 
element being collapsed on that scale. This corresponds to a dif- 
ferential mass df mass = p{S\R„) df mi = pc[Ri\ d f\o\- Since the total 
mass associated with the mass function is Mtot d/V (M) / AM, we have 
the predicted mass function or "first-crossing distribution": 

An _ pc{M) df 
dM ^ 



(21) 



M dM 

where d//dM is the differential fraction of trajectories that cross 5c 
between M and M + dM. 

This formalism has several advantages. It provides an exact 
solution that also allows us to rigorously calculate the normaliza- 
tion and shape of the mass function. It also allows us to explicitly 
resolve the "cloud-in-cloud" problem, i.e. to address the situation 
where a trajectory crosses the barrier multiple times. Figure [T| 
plots the resulting mass function (for a few choices of parameters, 
which just determine the normalization of the mass function and 
will be discussed below). We also compare the mass function if we 
were to ignore the "cloud in cloud" problem - i.e. where we treat 
every crossing above pc on a smoothing scale Rasa separate cloud. 
At the highest masses, the difference is small - this is because the 
variance is small and 5c is large, so the probability of being inside 
a "yet larger" cloud vanishes. However, at lower masses, the dif- 
ference rapidly becomes quite large (order-of-magnitude) - much 
larger than the factor = 2 of the Press-Schechter mass function. 
This owes to the complicated behavior of 5c, which increases again 
on small scales. Failure to properly account for the cloud-in-cloud 
problem and moving barrier will clearly lead to large inaccuracies. 

3.2 Key Behaviors 

If the barrier 5c were constant, the mass function of collapsed ob- 
jects would then simply follow the Press-Schechter formula; 

exp I - ^) (22) 



dnps po 2 5c \ diner | / i/" 
dM ^ M2 V ^ CT 1 dlnM V Y 
where v = 5c/(j{M) is the collapse threshold in units of the stan- 
dard deviation (cj(M)) of the smoothed density field on the scale R 
corresponding to M{R). 

However, the barrier here is not constant (it depends on R). 
A reasonable approximation to the first-crossing distribution, how- 
ever, is given by 



dn 
dM 



B = 



M^y TT a 

ln(pc.,niin/Po) 

ln(pc/po) 



1 dlno- 1 / 




IdlnMh^'H 


4) 



M < M{pc.„ 

M>M{pc,n 



(23) 



(24) 



crossings are not missed - i.e. so that the probability of crossing the barrier 
in a given step is small, A5 Sc{R). 



where pc. min = MW(pc[M]) is the critical density at the most- 
unstable scale. This is motivated by the exact solution for the first- 
crossing distribution for a linear barrier with Sc — Si+ a' /2, but 
with B — 5i held constant below M(pc)[^ Because the deviation 
from a constant barrier is only logarithmic, these formulae do not 
differ too severely, and we can gain considerable insight from their 
functional forms. 

Consider the behavior of both 5c and i/, which define three 
primary regimes. On scales above the sonic length but below ~ h, 
most of the dynamic range for GMCs, oc oc R"-^ (for power- 
law turbulent cascades) is large, so cr is a very weak function of R 
(most of the contribution comes from the largest scales, since p — 
1 > 0) while Pc decreases with R oc R''^^ so Sc — Inpc/ po ^ —{3 — 
p) InR. Therefore u (x 5c x -(3 - p) InR oc -((3 - p)/p) InM 
is a (logarithmically) decreasing function of mass. So we expect 
an approximately power-law mass function dn/dM oc M° with 
slope a ~ —2. This implies similar mass per logarithmic interval 
in mass and simply follows from gravity - which is self-similar 
- being the dominant force (since the turbulence is super-sonic). 
To the extent that the slope deviates from —2, it is because the 
barrier i/ gets larger towards lower M. From the above equation, 
M-^ exp ("i^V2) oc exp [- (3 - pf [In {M/Mo)f /2pV^] oc 
M" with 



aRi-2 + ln(Mo/M) (3- 
^ -2 + 0.1 log (Mo/M) 



pf/lp^a^ 



(28) 
(29) 



(where Mo is approximately the location of the mass function 
"break"; formally {A-K/Vjpoiv' ^ 10'' Mq for MW-like systems). 
In other words, we expect a slope a which is shallower than —2 by 
a small logarithmic correction, a ~ 1.7 — 1.9, as observed. 

At very small scales we approach the sonic length, M \; 
the growth in <7{R) becomes vanishingly small (~ \/3M/A oc 
^(p-i)/2-) ^yjjjig continues to increase logarithmically as before. 
The mass function must therefore flatten or turn over, with a rapidly 
decreasing mass in clouds below the sonic length (although the ab- 
solute number may still rise weakly). 

At large scales above ~ h, a{R) decreases rapidly with 
increasing R - the contribution from large scales goes as ~ 
^ln{l + {3/4)vj/K^R^) oc R-*+'' as -> oo, while now also 
increases oc R (so 5c oc InR), so the mass function is exponen- 
tially cut off as oc exp (— cM'"'''^'*). We caution that at the largest 
size/mass scales, global gradients in galaxy properties - which are 
currently neglected in our derivation of the collapse criterion - may 
become significant. However, the number of clouds in this limit is 
small. 



The fitting function from |Sheth & Tormen|j2002) : 

— dM = f(S)dS = \T(S)\exp[-5c(sY /2S] 



T{S) = J2 



(-5)" d"5c{S) 
~~n\ dS" 



(25) 



(26) 



gives a similar answer, but it is less straightforward to interpret. An approx- 
imate solution for the case neglecting the cloud-in-cloud problem is given 
by 

,,2 , 

(27) 



I do- 

- u 1 1 exp 

I dlnM" 



dn pc 3 ri dlnpr i 
dM ~ M2 y/2i(j y dlnM I 

which can be derived (up to a normalization) from differentiating Equa- 
tionfTsl 
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3.3 Comparison with Observations 

Figure [T| plots the predicted mass function: we show the exact so- 
lution, both excluding and including "clouds in clouds," and the 
approximations in Equation [23] & |27] For our "standard" model, 
we will assume the disk is marginally stable {Qo{h) = I), and that 
the turbulence, being supersonic and rapidly cooling, should have 
p ~2 (see the discussion in §[TJ. Motivated by observations, we 
normalize the turbulent spectrum by assuming a Mach number on 
large scales TVf;, ~ 30 (though we will show this exact choice has 
very weak effects, provided A4i, ^ 1). With these choices, the 
model is completely fixed in dimensionless terms. To predict an 
absolute number and mass scale of the mass function, we require 
some normalization for the galaxy properties: some measure of the 
local gas properties (mean density, velocity dispersion, surface den- 
sity, etc, to set the mass and spatial scales) and total galaxy mass or 
size (to know the gas mass available). Because of our assumption 
of marginal stability, many of these properties are implicitly related 
- we need only specify e.g. a total disk mass, gas fraction, and spa- 
tial size. Or equivalently, a mean density, velocity dispersion, and 
total mass. 

Taking typical observed values for the total gas mass, mean 
density, and velocity dispersion in the Milky Way, we plot the 
resulting predicted GMC mass function and compare to that ob- 
served. Because we are considering the total gas mass of the inner 
MW, we need to compare with a GMC mass function corrected to 
the same effective volume - we therefore compare with the values 
in [Williams & McKee| ( [1997] l (who attempt to construct a "galaxy- 
wide" GMC mass function for the same total volume). We then 
repeat the experiment with the average properties observed in the 
LMC and M33, and compare with the mass function compilations 
in |Rosolowsky| ( [2005^ ; |Fukui et al.H2008^ , corrected to the appro- 
priate survey area. 

In each case, the predicted mass functions agree remarkably 
well with the observations. We emphasize that although the ob- 
served densities and masses enter into the normalization of the mass 
function, the shape, which agrees extremely well, is entirely an a 
priori prediction. Moreover, the assumed densities do not entirely 
determine the normalization - because the barrier and variance are 
finite at all radii, the models here specifically predict that not all 
mass is in bound units. In fact, only ~ 20 — 30% of the total mass 
is predicted to be in such units - for the MW, the total bound GMC 
mass is predicted to be ~ IO'Mq, in good agreement with that 
observed jWilliams & McKee|1997^ . Likewise, the details of our 
stability and collapse conditions determine where, relative to the 
Jeans mass, the "break" in the mass function occursQ 

We should caution that it is not entirely obvious that our pre- 
dicted mass function is the same as that observed. The mass func- 
tion here is well-defined because we restrict to self-gravitating ob- 
jects and resolve the cloud-in-cloud problem, knowing the three- 
dimensional field behavior (and assuming spherical collapse). In 
the observations, the methods used to distinguish substructures 



' The predicted high-mass cutoff in the GMC MF is steep, but there is 
some suggestion that the GMC MF terminates or truncates more sharply at 



the maximum cloud mass in some systems (e.g. the MW; see [Williams &| 
[McKee|1997) . As noted above, including the corrections from global gra- 
dients in galactic properties in our collapse condition may steepen the pre- 
dicted cutoff. However, since the distinctions appear over a nan'ow range in 
mass (factor < 2) where the expected number of clouds in is in the Poisson 
regime (and consistent with zero within 2 a), it is difficult to discriminate 
between different models. 



and the choice of how to average densities (in spherical or arbi- 
trarily shaped apertures) can make non-trivial differences to the 
mass function ( Pineda et al.|2009) l. This may be considerably im- 
proved by the use of more sophisticated observational techniques 
that attempt to statistically identify only self-gravitating structures 
(see e.g. |Rosolowsky et al.| 2008); preliminary comparison of these 
methods in hydrodynamic simulations and observations suggests 
that most of the identified GMCs are indeed self-gravitating struc- 
tures so the key characteristics of the GMC MFs in our comparison 
should be robust, although details of individual clouds may change 
significantly ( [Goodman et aL|2009b| i. 

3.4 Effects of Varying Assumptions 

Of course, it is important to check how sensitive the predicted mass 
functions are to the assumptions in our model. Figure|2]shows the 
results of varying these assumptions. We plot the mass function in 
dimensionless units (po — h — 1, with the absolute mass being an 
arbitrary normalization). 

If we assume Kolmogorov turbulence (p = 5/3 instead of p — 
2), the predicted mass function is nearly identical at intermediate 
and high masses, but flattens more rapidly at low masses, because 
the velocity drops more slowly at small scales so pc oc R''~^ rises 
more steeply. The difference agrees well with the scaling in Equa- 
tion 28 which predicts a faint end slope a ~ —2 + 0.3 log (Mq/M) 
for p = 5/3 instead of a -2 + 0. 1 log (Mo /M) for p = 2. 

If we vary the Mach number on large scales A4h (or, equiv- 
alently, the assumed sound speed or magnetic field strength), the 
differences are very small at almost all masses, because the assump- 
tion that the disk as a whole is marginally stable effectively scales 
out the absolute value of A4i,- What A4i, does determine is the (di- 



mensionless) scale of the sonic length (/fsoi 



-2/ (;.-!) 



be- 
low which the mass function will flatten. With lower A4i, — 10, 
this happens at higher masses - but still quite low in absolute terms 
(R ~ 0.01 A, or > 3 dex below the maximum GMC masses). 

As noted above, the exact manner in which the velocity power 
spectrum E{k) should flatten at large scales < 1 is uncertain. 
We therefore re-calculate the mass function ignoring such flattening 
entirely ~ i.e. assuming E{k) oc k"'' for all k. This makes the very 
high-mass end of the mass function slightly more shallow, but has a 
negligible effect at all other masses. Since the only difference will 
be in the regime where the number of clouds is ~ 1 (so subject 
to large Poisson fluctuations) it is difficult to constrain this from 
observations. 

Re-calculating our results with a different window function 
makes little difference. We test this with a Gaussian window func- 
tion (convenient as it remains a Gaussian in real and Fourier space). 
As discussed in [Zentnerj ( [2007^ , this makes the calculation more 
complex because we can no longer treat the Fourier-space trajec- 
tory as having uncorrelated steps; following [Bond et al.[ ( [199l) the 
first-crossing distribution is computed by numerically integrating 
a Langevin equation. However, we hold our mass definition fixed; 
with this choice, for fixed R„, the exact choice of window shape 
about kw ~ l/Rw introduces only small (~ 10%) corrections (we 
refer to the discussion therein and [Maggiore & Riotto[[2010a[ for 
more detailed discussion of the effects of different window func- 
tions). 

What if the density distribution is not a lognormal? It has been 
suggested, for example, that for systems which have significant 
gas pressure and whose equations of state are non-isothermal, or 
which have large magnetic fields, the density distribution may more 
closely resemble a power-law (see e.g. ,Passot & Vazquez-Semadeni[ 
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[T998l|Scalo et al.|1998[[Ballesteros-Paredes et al.|2011b) . This is 
certainly still treatable with the excursion set formalism: there has 
been considerable discussion in the literature regarding the halo 
mass function and bias with non-Gaussian primordial fluctuations 



(see|Matarrese et aL|2000") [Afshordi & Tolley|2008| [Maggiore &| 
|Riotto|2010b| and references therein). However most of these rig- 
orous approaches assume the non-Gaussianity is small and can be 
treated in perturbation theory. For large deviations from Gaussian- 
ity it is not trivial to construct a fully self-consistent theory. For 
example, if P(p|i?„ ) were locally power-law at each "step" in k- 
space in a random walk, the resulting P{p) evaluated on each scale 
would no longer be a power-law; some violation of locality would 
be required so the distribution could "self-correct." In any case, if 
we simply assume some pre-specified /'(p|i?„ ) at all scales, it is 
still straightforward to evaluate the first-crossing distribution. The 
distribution f{S) in (d//dM) dM = f{S) dS is given by the solution 
to the integro-differential equation: 

m--P'iS.\S)'^-r'^dS (30) 



d5 



dS 



P'{5\S)=P{5\S)- r dS'fiS')P{5-S,[S']\S-S') (31) 
Jo 

where dp{S) — P{S)d5. This is essentially just the collapsed mass 
given by P{5 > 5c\S), corrected by the probability that the collapse 
occurred on a larger scale (smaller 5), and can be solved numeri- 
cally for any specified P. 

Consider the following form for the density PDF: 



dp(p) 
dlnp 



oc exp(-7|ln[p/p]|) : 



{p/py 

{plpy 



P<P 
P>P 



(32) 



where p = (1 — 7 po- The exact functional form is arbitrary, of 
course, but convenient because it is a pure power-law symmet- 
ric in Inp, and has a well-defined variance: {(Inp)^) = 27^^ + 
(In [1-7-^])- We can therefore map this one-to-one to our as- 
sumed density power spectrum by assuming 7 = y{P), with 
27~^ + (ln[l ~j~-])- = a^{R) =5. Note that this gives 7 ~ 1 over 
much of the dynamic range of interest, quite similar to the best-fit 
distributions in the references above. At low and high masses, the 
predicted mass function is slightly more shallow than our standard 
model. At high masses this is because of the more extended power- 
law tail to high S; at low masses this is both an effect of more first- 
crossings at larger scales and a result of some of the mass being 
moved from the "core" of the distribution to those tails. However, 
the differences are quite small. This is because a lognormal (unlike 
a pure normal distribution) is very similar to a single power law 
over a wide dynamic range. Moreover, the collapsed mass fraction 
is not extremely small, so it is not sampling some extreme tail of the 
distribution. So, for the same variance S, deviations from lognormal 
behavior have only small effects. 



3.5 The Core Mass Function 

In this paper, we choose to focus on the mass function of GMCs and 
other large-scale structures in the ISM. Part of the reason for this 
is that we can focus on the first-crossing distribution (the largest 
scales on which structures are self-gravitating) and so have a well- 
defined mass function. Although there are certain similarities, this 
is not the same as the mass function of self-gravitating dense cores 
within GMCs, as calculated using qualitatively similar arguments 
in e.g. [Padoan & Nordlund| ((2002) and [Hennebelle & Chabrier] 
l |2008) . 



In principle, our model can be extended iteratively to smaller 
scales to investigate the mass function of cores and make a direct 
comparison with these previous predictions as well as observations, 
and in a companion paper jHopkins]|2012^ we attempt to do so. 
This is not trivial, however. The difficulty is that, because cores 
are substructures, the mass function definition (the resolution to 
the "cloud in cloud" problem) is somewhat ambiguous: we can- 
not simply isolate first-crossing. Even in simulations where the full 
three-dimensional properties are known, it is not trivial to find a 
unique mass function of such substructure in a turbulent medium 
(see e.g. |Ballesteros-Paredes et al. 2006, Anathpindika 201 1). The 
approach of Hennebelle & Chabrier ( 2008 1 is to treat this ambigu- 
ity as an effective normalization term (and to truncate the problem 
at larger scales - treating the properties of the "parent" GMC as 
assumed/given and restricting to much smaller spatial scales); as 



such their derivation is similar to the original Press & Schechter 



|1974 ) derivation as discussed in §[T] That in [Padoan & Nordlund 



[2002 more simply makes some general scaling arguments. But as 
we show in Fig.[T| this is not necessarily a good approximation. We 
therefore require some more detailed criteria to inform our defini- 
tion of cores, for example some estimate of the scales on which 
fragmentation below the core scale will not occur (defining the 
"last-crossing," as opposed to "first-crossing" distribution). This is 
a topic of considerable interest, but is outside the scope of our com- 
parisons here. 



4 SIZE-MASS & LINEWIDTH-SIZE RELATIONS 

We can also use our model to predict the scaling laws obeyed by 
GMCs "at collapse." 

The linewidth-size relation follows trivially from our assumed 
turbulent power spectrum. The exact relation is plotted in 

Figure [3]for power-law turbulent slopes of p = 5/3 and p — 2, with 
the normalization set by requiring a marginally stable disk with 
MW-like surface density. We can define the line width either as 
just the turbulent width, or the turbulent width plus the contribution 
from disk shear a-l{R) — v," + R^; the distinction is unimportant 
for typical observed scales, but shear is predicted to contribute sig- 
nificantly to the velocities of the largest GMCs when R > h. We 
compare with obser vations compi l ed fro m th e MW and ot her local 



Bolatto et al. 



1 2008| and |Heyer et aL] ([2009i)p°] 



group galaxies from 

In the regime above the sonic length and below the scale 
height, this is just a simple power law with (T,,(if) oc i.e. 
~ 0.5 for p = 2 or ~ 0.33 for p = 5/3. This is essentially an as- 
sumption of our model (although it follows from basic turbulent 
conditions); more interesting is that the normalization can be pre- 
dicted from the assumption of marginal stability (Q « 1), giving 



a,.(.R) ^0.4kms" 



(Sdisk) 

10Mqpc-2 



XT' 



(33) 



This agrees well with the observed relation. In the full solution, 
because of the change in dimensionality above the scale h, the rela- 
tionship flattens if we consider only turbulent velocities; it becomes 
steeper, however, with the inclusion of the shear term. 

This model also specifically predicts a residual dependence 



Heyer et al. |j2009t caution that more detailed studies in nearby 
Goldsmith et al.|2008) suggest their LTE masses may be low by 



Because 
clouds (e.g. 

a factor ~ 2 — 3 at intermediate column densities, we plot the results for the 
"high density" cuts in cloud area defined therein (the "A2" sample) within 
which the LTE approximation should be valid. 
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Figure 3. Predicted GMC linewidth-size relation. Different lines corre- 
spond to different model assumptions: specifically we vary the turbulent 
spectral index (p), the absolute normalization of the system (amounting to 
the velocity dispersion cl at scale = h), and whether or not we include 
disk shear in the "velocity" a^.. Note that in the models here, ai, is not 
freely varied, it is predicted from the global parameters of the system via 
our marginal stability assumption. The velocity is the one-dimensional 
linewidth (using cr^ = Cj + v,-) for each cloud at the time of collapse, R 
is the three-dimensional collapse radius. On scales below ~ h, the Monte 
Carlo results are approximately a power-law with slope ct,, oc Ff ^ (Eq. |33| . 
We compare observations of clouds in the MW and local galaxies, compiled 
in |Bolatto et al.| j2008 circles) and Heyer et al. (2009 squares), appropri- 
ately corrected to the same quantities. The agreement is good - even for 
p = 5/3, for which large-scale effects make the relation slightly steeper 
than the naive expectation ay oc /f(''~')/2; moreover the marginal stability 
assumption predicts the normalization accurately. We also compare indi- 
vidual high-redshift molecular "clumps" in extremely gas-rich, rapidly star- 
forming lensed galaxies in |Swinbank et al.]^201 1| crosses with error bars), 
which form in much more dense disks (much larger Sdisk); these lie well 
above the extrapolation of the relation for MW-like properties. However, if 
we compare the predictions for a model with the observed cr;, ^ lOOkms"' 
of their host disks, the agreement is good. Clouds in the MW center, which 
has intermediate Sjisk between these extremes, lie correspondingly between 
these curves (see |Okaetal.|2001> . 



in the normalization of the linewidth-size relation that scales as 
(Edisk)''^", where (Edisk) is the large-scale mean disk surface den- 
sity. We stress that this is not necessarily the same as a dependence 
on the local cloud Ecioud (over a wide dynamic range, in fact, Edisk, 
hence Ecioud, is similar). This is also, by definition, for bound ob- 
jects, not for un-bound overdensities on small scales. The predicted 
dependence is shown indirectly in Figure [3] and directly in Fig- 
where we compare with the observations compile d in | Heyer| 
^2009^ in local galaxies and in [Swinbank et al.| ( |201l| for 
massive star-forming molecular complexes in lensed, high-redshift 
galaxies. These sample extremely different environments, and are 
indeed offset in the linewidth-size relation. However, the magnitude 
of their offsets is in good agreement with that predicted here The 
galaxies in [Swinbank et alT]j2011^ have an average surface den- 
sity of ~ 10''Mqpc~'^, and a correspondingly very large measured 
ah ~ 100 kms^ (as expected for Qo{h) ~ 1); normalizing the pre- 



" If the predicted clouds perfectly followed M oc and a oc R'/-, they 
would collapse to a single point in this Figure. They do not, because of the 
changes below the sonic length and above ~ h. However, because the clouds 
are defined as self-gravitating, the models collapse to a line (with most of 
the clouds concentrated near the "typical" point for intermediate scales. 
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Figure 4. Size-mass relation of clouds, in the same style as Figure |3] 
The observations from |Bolatto et ar]j2008) us e th e virial mass estimator 



Mvir = 5 (jIR/G; those from Heyer et al. ' 2009| and 



Swinbank et al. 



j201lt 



are derived from the CO luminosity. The agreement with observations is 
good, and the scaling is an approximate power-law with slope R oc M"'' 
(approximately constant Scloud ^ IOOMq pc^^; Eq. |34| . Again, the high- 
redshift clumps lie off the "typical" local galaxy relations; however a model 
of a more dense disk with c/, Ki 100 agrees well. Similarly, MW-center 
clumps lie between the extremes shown ^Oka et al.|2001) . 
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Figure 5. Residuals from the linewidth-size relation for hound clouds as a 
function of disk/region surface density S, in the style of Figure]!] Because 
we define clouds as self-gravitating, the different predicted lines (from dif- 
ferent turbulent spectra) in Figs. |3|4| collapse to a single line in this plot. So 
we instead plot the predicted lines for an assumed global stability parameter 
2o ~ 0-5 — 2.0. Unbound clouds/overdensities will have higher a,, but are 
not the collapsed objects followed here. 



dieted linewidth-size relation for these properties, we expect an or- 
der of magnitude larger cr,, at fixed size. Clouds observed in the MW 
center ( |Oka et al.|2001^ , which has a higher mean surface density 
than the local neighborhood but generally lower than estimated for 
the high-redshift systems, lie neatly between the predicted curves 
for the local and high-redshift cases (a mean offset of ~ 3 — 5 rel- 
ative to the local clouds, corresponding to a factor of ~ 10 — 30 
higher Edisk, about what is expected for the observed exponential 
profile). Similar offsets are known in other local galaxies with high 
surface densities, such as mergers and starburst galaxies ( |Wilson| 
|et al.|2003[|Rosolowsk y & Blitz 2005). 



As discussed in Hopkins et al.^ lj20T2J, a dependence of ex- 
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actly this sort is seen in highi-resolution hydrodynamic simula- 
tions as well. In the observations, this normalization dependence 
has sometimes been interpreted as a consequence of magnetic sup- 
port or confining external pressure (see the discussi on in [Blitz &| 
|Rosolowsky|2006| [Bolatto et al.||2008i [Heyer et al.||2009| , but in 
this context magnetic fields and pressure confinement are not ex- 
plicitly present - such a scaling is a much more broad consequence 
of the simple Jeans requirements for collapse in any marginally sta- 
ble environment. 

The size-mass relation follows from the critical density pc de- 
rived in § |2.1[ by simply inverting Eqn[T4] We plot the exact predic- 
tion in Figure|4] In the regime above the sonic length but below the 
disk scale-height, recall that a power-law turbulent cascade gives 
the simple condition pc = v,{k)- / {4tvG) oc R''~^, so oc M'^'', 
i.e. R oc M^^^ for p ~ 2, very similar to the observed power-law 
scaling. The normalization also follows - for MW-like global con- 
ditions 



10 



^doud(^>4omc,^</!) 



1.40-0 4' pc I 



V 300 Ms,,,,/ 



0.5 



(34) 



where aoA is the normalization of the turbulent velocities v, — 
<JoA X 0.4kms^' (R/pc)'^^. This corresponds to an approximately 
constant cloud surface density in agreement with Larson's laws: in 
projection Ecioud ~ IOOMq pc~^ at the time of collapse. Note that 
re-calculating this for p = 5/3 only changes the slope from 0.5 to 
0.6, well within the observational uncertainty. This will also alter 
behavior at the highest masses, but this is not significant until well 
above the mass function break. There does however appear to be 
tentative evidence for such a transition in the observations shown 
in Figure [4] As expected from the behavior of the linewidth-size 
relation, clouds in high density environments - which will have a 
higher cro.4 in Eqn.[34]above - are offset to lower R at fixed Mdoud; 
we show the same model prediction for the high-redshift systems 
in good agreement with the observations. Once again, MW center 
clouds and other local systems in environments with higher densi- 
ties are similarly offset. 

As discussed in § |3.5[ fully extending the models here to the 
scales of dense cores is beyond the scope of this paper. However, 
we expect these cores, if self-gravitating, to obey the scaling in Fig- 
ure |5] This means that if they form inside of high-density GMCs, 
we can (approximately) think of the "parent" GMC surface den- 
sity as similar to the background (Edisk) term in Eqn. 33 and 
might expect them to have higher dispersions at fixed sizes. This 
has been seen in suggested from observations jBallesteros-Paredes] 
|et al.|2011a^ , as part of a quasi-hierarchical gravitational collapse, 
similar to the predictions here. Of course, some regions can have 
much higher cr,, at fixed R and be simply not self-gravitating; these 
will not lie on the relation in Figure |5] (they will be offset to 
higher ct,,//?'''^). This may, in turn, give rise to a dependence of 
the linewidth-size relation on the tracers and extinction threshold 
adopted, as observed ( [Goodman et al.|199F{|Lombardi et al.,2010^ 



5 SPATIAL CLUSTEMNG OF GMCS 

In analogy to dark matter halos, we can use the excursion set 
formalism to also determine the spatial clustering and correlation 
function strength of these bound sub-units. Following |Mo & White| 
( [1996^ , the excess abundance of collapsing objects (relative to the 
mean abundance) in a sphere of radius Ro with mean density So is 



Scon{Ri, Sc.i \Ro, So) 



AA(1|0) 
n{Mi)Vo 



1 



(35) 
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Figure 6. Predicted linear bias b - i.e. the amplitude of .spatial clustering - 
as a function of GMC mass (allowing for clouds inside of bound over densi- 
ties). We plot this for our standard model and model variations in Figure|2] 
in dimensionless units. Low-mass GMCs are weakly biased or anti-biased 
- they simply trace the dense gas. The highest-mass GMCs are strongly 
clustered - they preferentially trace global overdensities (e.g. spiral arms, 
galaxy nuclei, etc.). 



where n{Mi) is the average abundance of objects of mass Mi (from 
the mass function) and Af{ 1 10) is the number of collapsing objects 
in a region of radius Ro (variance So) with fixed overdensity So- 



5.1 Linear Bias 

If Sc were constant, A/'(1|0) can be determined analytically and is 
simply 



mm- 



Pc. 1 Vo Sc. I 



M 



exp 



{ScA-Sofl d5i 



1 V27r(5i-5o)3/2 'L 2(5i-5o)JdM, 

(36) 

( |Bondetal.|1991^ . In the regime where > ^ 1 , so Ao < A 1 , this 
simplifies to 



Scot 



So = h(M{)So 



(37) 



where b[M\ ) is defined as the linear bias of objects of mass M\ p] 
The barrier Sc here is not constant. However, for arbitrary 
5r(M), we can calculate A/'(1|0) exactly by repeating our Monte 
Carlo excursion from § |3.1[ but instead of beginning with initial 
conditions S = 0, J = for each walk, we begin at scale S = So 
with density S = So- The bias b{Mi) is then just the ratio of 5coii/5o 
for small Sq. 

Figure [6] plots the bias as a function of cloud mass. A couple 
of key properties are clear. At high masses above the exponential 
cutoff in the mass function, the bias increases rapidly. This is qual- 
itatively similar to what is seen for dark matter halos: because such 



'■^ The expression for bias here is different from that for dark matter halos 
by a linear offset of unity. That offset arises in the dark matter case because 
of the expansion of the Universe and subsequent mapping from "initial" 
(Lagrangian) coordinates to "observed" (Eulerian) coordinates. It does not 
appear here because the terms are all evaluated instantaneously (the expres- 
sion here is equivalent to the "initial time" expression for b in halos). 
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Figure 7. Comparison of the predicted correlation function ^(R) of bound 
gas objects/GMCs (lines) with the observed congelation functions of young 
star clusters. We show the predicted values for three different masses (essen- 
tially a normalization difference in the correlation function, corresponding 
to the bias b in Fig.|6j. For lower masses, b changes weakly, and for higher 
masses the MP is exponentially suppressed, so this covers the interesting 
range. We plot radii in units of the scale height h, in which the correlation 
function is dimensionless. Top: i;^^, the cross-correlation between bound 
objects and gas mass (defined in Eq.|40j. We compare the observed cross- 
correlation between young star clusters (which should trace the locations - 
regardless of how efficiently they form - of their "parent" GMCs) and CO 
gas maps measured in the Antennae by |Zhang et al.|{2001) . We compare 
the two youngest ^ 100 cluster samples (sampling two different regions 
and mass ranges), with ages < 5 Myr and ~ 3 — 16 Myr We do not know 
the masses of the progenitor GMCs, but they are likely to be in this range, 
since these are the most massive young star clusters in the galaxy (the more 
massive sample has the higher l^cm |). Despite this being a disturbed system, 
the agreement is reasonable. Bottom: ^mm, the auto-correlation function of 
bound objects. We again compare this measured for the young star clus- 
ters in the Antennae (orange circles). We also compare the young ~ 1000 
star cluster autocorrelation function in M5 1 (which is not disturbed), mea- 
sured by |Scheepmaker et al.H2009) for age intervals 2.5 — 10, 10 — 30, and 
30 — 300 Myr (red, cyan, green, respectively). Especially in the youngest 
samples, the agreement is good. We compare the same, measured directly 
for GMCs with mass ^ 2 Sj /i^ in M33 from |Engargiola et al., ( ^2003^ ; again 
the agreement is good. 



systems are exponentially rare, they will tend to be strongly bi- 
ased towards the few regions with substantial large-scale over den- 
sities. Physically, this corresponds to gas overdensities in the disk 
on scales larger than the scale-height h, i.e. a preferential concen- 
tration of the most massive GMCs in global instabilities such as 
spiral arms, bars, and ~kpc-scale massive star-forming complexes, 
rather than their being randomly distributed across the gas. At in- 
termediate masses below the mass function break, where most of 
the cloud mass lies, the bias is weak (order unity), so most of the 



mass in clouds simply traces most of the gas mass in general. We 
stress that this does not necessarily mean clouds are randomly dis- 
tributed over the disk as a whole; it means they are unbiased relative 
to the gas mass distribution. But at low masses, the bias again rises 
(weakly). This is related to the anti-hierarchical nature of cloud 
formation: the bias here is driven by clouds which form via frag- 
mentation from other clouds. 

We can approximate these exact results using our previous ap- 
proximate fitting functions for the mass function (Equation [23] & 
\27\ modified (as with the case of a linear barrier) so S — > 5i — 5o and 
u {5c — Sof / {Si — So)- Neglecting the clouds-in-clouds prob- 
lem (i.e. including those clouds), we obtain the approximate 



5cA{i + Vi dlnpi/dai) 



(38) 



Which in practice is a small (~ 10 — 20%) correction to Equa- 
tionlSVl If we exclude clouds-in-clouds, 



b{Mi 



1 r 2 Sc.n 

^i r - 



(39) 



(where B is defined in Equation 23 i. This is identical to Eqn. 37 
at high masses, but it allows for negative bias at low masses, if 
Sniin = ln(pr.mm/po) < 2 and Sc < o- (S^,[, - 1/2). Physically, the 
fact that Equation [38] is always positive means that the number of 
bound regions of mass Mi inside a large-scale overdensity always 
increases with 5o. However, for some values of Mi and So, increas- 
ing &Q more rapidly increases the probability that these regions are 
themselves inside a larger collapsed region. For a more detailed dis- 
cussion of the leading-order corrections when considering a mov- 
ing as opposed to constant barrier 5c, we refer to |Sheth et al.| ( |200"T| l. 

5.2 The Correlation Function: Theory 

Recall, the physical over-density is p/ po = exp {5 — a{R)^/2). The 
correlation function between collapsed objects of mass Mi and 
background mass, as a function of radius Ro, is defined by 



l+Ccm(^0,Ml 



^ {m\Ro)\p) 

n{Mi)Vopo 
= ((1 + 5coii) 1 exp {So - S[Ro]/2))r„ 

4M-.(^°-^°/^',(5o|5o)d5o 
n{Mi)Vo 



(40) 
(41) 
(42) 



where the integral is over all So < Sc{Ro), and q{So \ So) is a weight- 
ing factor defined in |Bond et al.| ( fT99Tj l as the probability that the 
overdensity at a random point, smoothed on a scale Ro, is So and 
does not exceed Sc{Ro) on any larger smoothing scale. 

Equation |42] can be evaluated numerically with the Monte 
Carlo solution for A/'(1|0) and q{5o\So). But at large 3> Ri (pro- 
vided 5o — >■ as ^ — >■ oo), it simplifies to just 

l+^cr.{Ro,Mi)^i+b{Mi)a-{Ro) {Ro»Ri) (43) 
= l+b{Mi)UJ.Ro) (44) 

This can be shown for any first-crossing distribution by first tak- 
ing q — > p{Sq 1 5o) since the probability of collapse on larger 
scales is negligible, and then noting exp{5o — So/2) p{5o\So) ~ 

Note that the equations here are modified from those used in the cos- 
mological case because we use 5 to represent the logarithmic (not linear) 
density field. However for small Sq they are identical, which is why we 
recover similar scalings for the bias and correlation functions. 
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1 / \/2vr5o exp {—{So — 5o)"/25o), which becomes a deha function 
centered at So — So as 5o — )■ 0. 

The auto-correlation function of the mass ^mm is given by 
1 + Cmm = (p^) /pI = exp (5o), so ^mm = exp (5o) - 1 ^ So = cr^ {Rq) 
at large Rq is just the variance in the mass field. So the collapsed 
object-mass correlation function on large scales is then just the bias 
times the mass autocorrelation function. It is straightforward to ver- 
ify that the auto-correlation function of collapsed objects is just 
given by 



Ccc~fe (Mi)e.m {Ro:>Ri) 



(45) 



The correlation functions discussed above are the three- 
dimensional correlation functions. However, with rare exceptions, 
it is in general much easier to determine the projected correlation 
function ^2d{Rp), defined so the probability of finding another ob- 
ject in a 2-d annulus d~r around a given object is (dN/dA) (1 + 
(_2d)d^r. This is straightforward to calculate 



r^no{z)^3,{VW+?)'i^ 



(46) 



where z is the line-of-sight direction and no{z) = n{M) is the av- 
erage abundance. For the typical case of an approximately face- 
on disk with the exponential vertical profile we have adopted, 
wo(z) oc exp(— z/A); however, accounting for this, we should also 
slightly modify our calculation of ^3^, integrating over po at all 
central positions with A/'(l |0|x) a function of po(x) (since our 
derivation up to this point implicitly assumed a homogenous back- 
ground). In either case, at large radii this is just w,, oc {R/h) ^jd. 



5.3 Observed GMC & Star Cluster Correlation Functions 

In Figure [7] we compare the predicted (two-dimensional) correla- 
tion functions to observations. Unfortunately, at present there are 
no published observations of the GMC-GMC correlation function. 
However, various groups have measured the correlation functions 
of young, massive star clusters in nearby systems. Statistically, the 
positions of such star clusters should trace those of their "parent" 
GMCs (with greater fidelity as we consider younger star clusters). 
And although clusters will disperse or be destroyed with time, the 
correlation function should not be affected so long as this "infant 
mortality" is not strongly position-dependent (though that is uncer- 
tain, if it depends on e.g. tidal fields). This also has the advantage 
that star clusters can be much longer-lived than GMCs, so allow 
better statistics. The major uncertainty is that, without knowing the 
(uncertain) star formation efficiency, the exact mass of the progen- 
itor GMCs is undetermined. However, since the observed systems 
sample the brightest clusters, we can safely assume that their pro- 
genitors were the most massive GMCs (and since the mass func- 
tion cuts off exponentially, should reflect masses ~ 1 — 10 times 
the "break" in the mass function). 

|Schee pmaker et al.|i2009 1 measure the star cluster-star cluster 
auto-correlation function (which we should compare to the GMC 
autocorrelation ^cc) in M51 for the brightest ~ 1000 star clusters, 
in three age intervals (2.5 - 10, 10 - 30, and 30 - 300 Myr). The 
cluster masses range from 10'''~^ Mq, which for a few percent 
star formation efficiency indeed corresponds to the most massive 
GMCs. The mass scale only affects the bias (normalization) - it is 
more important to compare the shape of ^cc - this is invariant in 
units of R/h. With a large number of clusters and a nearly face-on 
projection, this is the most robust probe over large dynamic range. 
[Zhang et al.| ( |200T] l measure in the Antennae the star cluster-star 



cluster autocorrelation function and the star cluster-gas cross cor- 
relation function (tracing the gas in CO maps, which - since the 
system is quite dense - account for most of the gas mass). Here the 
geometry is obviously much more complex so the results should be 
interpreted with additional caution, but the authors do attempt to 
account for the global structure of the system, and separately mea- 
sure the correlation functions in different regions. We specifically 
consider their youngest cluster samples (R and B 1 ), with the bright- 
est ~ 100 and ~ 1000 objects at ages < 5Myr and ~ 3 — 16Myr, 
respectively (masses ~ 10'*~*Mq). Finally, we attempt to follow 
the procedure in |Scheepmaker et al.| (2009) to construct the auto- 
correlation function for GMCs in M33, using the catalogue in |En-| 
Igargiola et al.| ( (2003] l, which is both face-on and has a well-defined 
survey area and completeness limit. Since we cannot properly ac- 
count for survey edge effects or the global density profile, we sim- 
ply truncate the correlation function at half the radius inside of 
which 75% of the identified GMCs are found. Here, we can de- 
termine the mean mass in the distribution, which is approximately 
~ ITihh' estimated using the parameters from Figure [I]- this is 
almost exactly the value in the model which gives the best-fit pre- 
dicted normalization of S,{R). Given the uncertainties in both obser- 
vations and the cluster-GMC mapping, the agreement is striking. 



6 THE DISTRIBUTION OF UNDERDENSE BUBBLES 

Just as we used the excursion set formalism to predict the mass 
function of clouds by identifying objects above a critical over- 
density Sc, we can also use it to predict the abundance of under- 
dense regions ("bubbles") by identifying regions below a critical 
under-density Si,. We will follow [Sheth & van de Weygaert| ( (2004| , 
who apply this formalism to the dark matter halo context to study 
the distribution of voids. 

Generally, the procedure is the same, but considering the 
mass/radii below Sb instead of above Sc. However, some additional 
complications arise. First, unlike the case of collapsing objects 
where the counting of "clouds in clouds" was potentially valid, 
here we should clearly count "voids in voids" as simply part of 
the larger, parent void/bubble. So we again need to specify to the 
first crossing distribution (the distribution of the largest radii on 
which trajectories cross Si,). Second, we must also ensure that the 
void/bubble region is not itself contained inside of a collapsing re- 
gion (i.e. that S < Sc on all scales above the Sb crossing), since that 
would "overwhelm" or "squeeze' ' the bubble[3 Third, and most 
critical for our purposes, a "void" or "bubble" is not obviously well- 
defined in this context. Because there is no linear expansion here, 
we cannot derive the equivalent of the shell crossing criterion used 
for dark matter halo voids, and there is no obvious threshold which 
is physically as robust as the self-gravity criterion for collapse. We 
will return to this question and consider different plausible, but ul- 
timately somewhat arbitrary choices of under-density criterion. 

If the "bubble" barrier Sb and the collapse barrier (which must 
be avoided on scales above the bubble) Sc were constant, then |Sheth| 
|& van de Weygaert| ( |2004| l show that the first crossing distribution 
can be analytically re-derived subject to these boundary conditions. 



The details of the criterion for this can be subtle and more complex than 
simply being in a collapsing region, since smaller overdensities can also 



"squeeze" voids. This is discussed in detail in |Sheth & van de Weygaert] 
l2004i. However, because we do not need to map here between initial and 
final overdensities, many of these ambiguities are avoided. 
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Figure 8. Predicted differential volume fraction in underdense "bubbles" 
as a function of bubble radius R. For illustrative purposes we assume 
h = 200 pc but the scale Sbubblc scales oc Top: Bubbles defined as a pro- 
portional under-density p < po/100- Bottom: Bubbles defined as regions 
where the post-shock cooling time at velocities ~ v, (R) exceeds the free- 
fall time ~ 1 / \/Gp. Because determining the cooling time requires absolute 
units, we normalize the model by assuming S = IOMqpc^^. The exact 
solution (black solid lines) is compared to the approximate analytic solu- 
tion (red dashed) from Eq.[^ A broad distribution of underdense regions 
should be present simply from turbulent velocity divergences, which can 
have sizes ~ h and contain a large fraction of the disk mass. However, only 
a small fraction will "self-heat" to temperatures where they cannot cool - 
"hot" bubbles require energy input from some source (e.g. stellar feedback). 



to give the fraction of trajectories in bubbles per logarithmic inter- 
val dlnui, 



Vbfb{vt) = ^ ^"^2^ sin(;j7rr>) exp - " ^ 



V = 



Vb 



(47) 



(48) 



'■'-5(.R)'/2 

Recalling that we are sampling the Eulerian space, we can then 
trivially translate this to the number density of bubbles per unit 
radius or mass, e.g. 

An _\ df 
^ Vb d\nR 



1 dlnz/f, 

■ Vhjb(yh) — 



(49) 



dln.R Vb d\nR ' ■ ' Vb d\xs,R 
where 14 is the effective volume of the bubble. 

Again, we stress that the barrier is not constant, so we do not 
know that this will be an accurate approximation. More rigorously, 
it is straightforward to derive the same first-crossing distribution 
using the Monte Carlo approach in § |3.1| We follow the identical 
procedure, but simply record the first crossing of &b(R) for those 
trajectories that cross < <5i(i?) and have not crossed &c{R) at 
any larger scale. 
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Figure 9. Comparison of predicted and observed hole/bubble radii. We plot 
the predicted cumulative number of bubbles as a function of bubble size 
for our standard model, in dimensionless units (bubble size in units of A). 
Here, we assume a simple order-of-magnitude proportional bubble under- 
density p < po/10- For typical galaxy properties, this also corresponds to 
densities at which the diffuse galactic UV background will fully ionize the 
bubble. This allows us to plot all observed systems on the same Figure. We 
compare the observed HI hole radius functions from radius functions from 
the S MC jStaveley-Smith et al.|1997> , Holmberg II jPuche et al.(T992j , and 
M31 ^Brinks &^^^a[l986^ , and normalize them with the observed Mgas, 
po, h from the same sources. The agreement is good - most, if not all, of the 
HI "holes" are a natural consequence of turbulent density fluctuations and 
require no input energy source to "clear them out." 



The results of this exact calculation, and the analytic approx- 
imation from Equation |47] are shown in Figure |2] for two differ- 
ent choices of ib- First, we consider a simple under-density cri- 
terion: here p], < po/100. There is a very broad distribution of 
bubbles which satisfy this criterion: it includes several tens of per- 
cent of the total mass. The characteristic spatial "bubble scale" is 
at a factor of ~ 0.1 A, which (for the definitions used here) corre- 
sponds very closely to the scale at which the local contributions 
to density fluctuations (AS) are maximized. A large population 
of such fluctuations must arise for a density distribution similar 
to Equation |6] because the distribution is lognormal, the median 
density is ln(p„,cd/po) = —a'/l; i.e. for a ~ 1.3 dex fluctuations, 
Pracd « 0.01 po so of order half the volume should be in underdense 
regions. For any fixed (fractional) density threshold pb/po, the be- 
havior is qualitatively similar, but shifts systematically to smaller 
scales R and smaller normalization (the total mass in such regions 
scaling as ~ exp (— t'j /2) as pb / po decreases. 

There is nothing physically "special" about such regions - 
they are simply the low-p portion of the density PDF. A more mean- 
ingful threshold might be to define "bubbles" as regions where the 
cooling time becomes longer than the free-fall time. The isothermal 
temperature Cj is however quite low, so this will not be satisfied 
unless the temperature suddenly increases; for this, consider the 
shocks occurring in the turbulent medium at v., ~ vv (R) . Knowing 
E{k), we can estimate the distribution of post-shock temperatures 
and densities for a random Lagrangian parcel, and compare the re- 
sulting cooling time to the free-fall time fff ~ 0.54/ y^G p. Since we 
are interested in the regime where the cooling time will be long, 
we can simplify the problem by assuming a strong adiabatic shock 
and that thermal Brehmsstrahlung dominates the cooling. In this 
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regime, fcooi = nksT /An' < tff at densities 



(50) 



If we normalize our model to MW-like conditions by assuming |1974[ [Williams & McKee|[T997l|Evans|1999nEvans et al.,2009 
ag{h) ~ lOkms^' and no = po/fJ-mp ~ lcm~^, then this defines 
Si,. The resulting distribution of bubbles is shown in Figure[8] Qual- 
itatively, the shape of the distribution is similar - it truncates more 
rapidly at low R because the decrease of turbulent velocities (v, (R)) 
with decreasing R means that the barrier becomes more difficult to 
cross. The normalization is also significantly lower, corresponding 
to the lower absolute densities {pb/po ~ 10^*) needed near scales 
~ /i to reach this "hot gas" threshold. 

In both cases, the analytic approximation of Equation [47] 
works well for the largest voids (albeit with a factor ~ 1 — 1.5 nor- 
malization offset), but is systematically offset for low-mass voids. 
This is directly a consequence of the moving barriers 5c and Si,. 

In Figure|9] we compare the predicted size function of bubbles 
(in dimensionless units) to observations of HI "holes." We compile 



the observed HI hole size distributions in the SMC iStaveley-Smith 



eraTl[T997l, Holmberg II ( [Puche et al.][l992 l, and M3I ( [Brinks 



& Bajaja 1986), and scale the normalization of each according to 
the observed global galaxy properties measured at the radii enclos- 
ing half the "holes." Observations of the LM C ([Kim et al.|1999l), 
IC2574 ( [Walter & Bniiks[[T999) , and M33 ( [Deul & den Hartog[ 
[I990| l give similar results. 

There is no well-defined criterion for selection of "holes" and 
the density contrasts involved are typically modest, so we sim- 
ply compare with the prediction for a constant density contrast 
Pb < po/10. This is approximately consistent with the direct es- 
timates of the interior bubble densities/density contrast, and also 
(for the global properties of these systems) corresponds to densi- 
ties where even the largest (few hundred pc) holes would become 
fully ionized either from the diffuse galactic background or a sin- 
gle 0/B star inside the "hole." The agreement is good - if anything, 
the model predicts more small "holes," but this may be a question 
of observational selection/completeness (or a deficit of sources to 
ionize them). The characteristic hole size is predicted to scale with 
h (the characteristic radius), giving larger holes in thicker galaxies 
- a well-observed phenomenon (see |Oey & Clarke [1997 [[Walter &[ 
[Brinks|I999[ and references therein). 



7 CONSTRUCTION OF GMC "MERGER TREES" FROM 
TfflS FORMALISM 

7.1 General Considerations 

One of the most powerful applications of the excursion set ap- 
proach in galaxy formation is the construction of the extended 
Press-Schecter "merger trees," conditional mass functions, and for- 
mation histories for dark matter halo populations, which form the 
foundation of semi-analytic models. This provides a means to sta- 
tistically link populations in time and self-consistently model their 
evolution, with whatever additional physics are desired. We now 
show that the same "merger tree" approach can be applied here, to 
derive the time evolution of the systems we have thus far consid- 
ered static. 

Before we describe the mechanics of constructing these trees, 
there are a couple of important physical distinctions that will neces- 
sitate a somewhat different treatment from the typical methodology 
in the dark matter halo EPS formalism. 



First, unlike with dark matter halos, there is no reason to be- 
lieve that bound clouds are "conserved" (modulo their mergers into 
more massive clouds). In fact we expect from observations that 
they only live a short time, then are disrupted ((Zuckerman & Evans] 

?l |Evans et al.,2009^ 

So it makes no sense to begin from a present population of clouds 
and work backwards in time to construct the tree (as is typically 
done for halo merger trees). Instead, we need to forward-model the 
time evolution, to allow whatever model physics the user desires to 
determine whether or not such clouds survive. 

Second, we cannot assume that all the mass is in collapsing 
objects. We must therefore track un-coUapsed elements as well, al- 
lowing for their possible collapse at later times. 

Third, density fluctuations in a turbulent medium clearly do 
not evolve according to simple linear growth, in the manner of 
cosmological density perturbations. How, then, can we link a fluc- 
tuation at any one time to that at another time? To do this, we 
will assume that the turbulence is globally steady-state: i.e. that 
- excepting the behavior of collapsing regions - the turbulent ve- 
locity cascade is (statistically) maintained and, as a result, so is 
the global density PDF. We stress that we are not attempting to 
model how the turbulence is maintained. In this regime, the density 
PDF for independent modes on different scales obeys a generalized 
Fokker-Planck equation, with a diffusion term giving the effectively 
"random walk" behavior of each Lagrangian density parcel (from 
small-scale encounters/shocks/accelerations) and a drift term cor- 
responding to damping/relaxation (from viscosity, pressure forces, 
mixing, and the energetic cost associated with large velocity devia- 
tions). Under these conditions, if we know the stationary behavior 
of the PDF for some variable x is a Gaussian distribution with stan- 
dard deviation a, and zero mean, then the probability distribution 
to find the system with value x at time t given an initial distribution 
with (delta-function) value xo at time /o = f — A ? is 



p{x, t) dx - 



1 



\/27ro-| 



exp 



(x — xoY 
2 a? 



(51) 



= Ox v/l-exp(-2[f-fo]/r) a,. ^/lEtJ^ (52) 
= xo exp (-[/ - ?o/r]) xo (1 - Ar/r) (53) 

where the latter equalities are the series expansion for Af/r< l[3 
The timescale r here is the timescale on which the variance 
of x(f) with respect to xo grows, normalized by the steady-state 
variance cr., , i.e. the timescale of "mixing" in the distribution. More 
formally the amplitude of the correlation between values in time de- 
clines with exponential timescale r. In supersonic turbulence, this 
is simply the crossing time 



(54) 



Where ~ I is a constant which can be calibrated from numerical 
simulations ( [Pan & Scannapieco,^2010, find 77 ~ 0.90 — 1.05 over 
the range ~ 1.2-10). 



In addition to being convenient later, these series expansions have the 
properties that for small timesteps, they represent the only form that the 
evolution of p{x, t) can take if we require that it conserve Cv in ensemble 
average and conserve the growth in variance between x and xq independent 
of the choice of integration stepsize. 
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1.1 The Mechanism of Tree Construction 

With these points resolved, it is straightforward to generalize our 
approach to construct a time-dependent "fragmentation tree." We 
outline the methodology below. 

(0) Define the variance 5 = <J^{R) and collapse threshold 
5c.{R) either directly or from the turbulent power spectrum E{k). 

(1) Begin by constructing the initial conditions. Consider a 
Monte Carlo ensemble of "trajectories," as in § |3.1| Each trajectory 
5{R) is defined by the values ASj on each scale Rj — >■ Rj — AR. 
We are free to choose whatever values of ASj define an appropri- 
ate initial condition. For example, we can assume that the medium 
has a density PDF corresponding to "fully developed" turbulence 
and generate the ASj exactly in § |3.1| Or we can begin with a per- 
fectly smooth medium, setting all ASj = 0, and treat all structures 
self-consistently as they develop. Critically, save the full trajecto- 
ries ASj (full 5{R)) for each member of the Monte Carlo popula- 
tion, including those for which the region is "uncollapsed" (S{R) 
never crosses Sc{R)). 

(2) Evolve the system forward by one time step At. For a 
Fourier-space tophat window, we can evolve the system by perturb- 
ing each ASj independently according to Equation |5T] (obtaining a 
new, perturbed trajectory S(R, t + At)). The probability distribu- 



tion for the perturbed A(5;(i?, / + A?) is given by EquationSl with 
the appropriate substitutions: 

Ap{ASj[t + At]) _ 



d(A(5,[/ + Af]) 
1 



(55) 



exp 



{ASj[t + At} - ASj[t}./T^Y 



a/Ztti/'AS 

V) = 1 -exp(-2Af/r) 

r^R/{v,{R)y 
This is equivalent to taking 

A5j{t + At) = ASj(t) exp(-A?/r) 



2^ AS 



(56) 
(57) 



(58) 
(59) 
(60) 



+ n ^AS ( 1 - exp (-2 At/r) ) 

^ ASj(t) (1 - Af/r)+7^v/2A5Af/r 

where 7?. is a Gaussian random number with unity variance. This is 
done for all ASj in the trajectory, giving a new trajectory 

Rj>R 

5{R,t + At)=J2^Sj{Rj,t + At) (61) 

which can now be evaluated. 

(3) After each timestep, evaluate all trajectories 5{R) in the 
Monte Carlo ensemble. If the trajectory did not cross 5c(R) at any R 
in the previous time steps - i.e. it represented an uncollapsed region 
- then either it will remain uncollapsed (5{R) < Sc{R) at all R) in 
the new time step, or it will now cross the barrier at some R^. The 
largest such Rc corresponds to the collapse scale, defining a new 
self-gravitating object with mass M = 47r/3 Pc{Rc)RI. Physically, 
this event corresponds to the random density fluctuations from e.g. 
shocks and other processes pushing this previously "diffuse" gas 
parcel to densities at which it becomes self-gravitating, and col- 
lapses. The trajectory should still be saved, but the mass is now 
in a self-gravitating object, and the first-crossing scale on which it 
became self-gravitating should be recorded. 

If the trajectory already crossed the barrier at some R^, then 
there are two possibilities. If the trajectory no longer crosses the 
barrier (or crosses at some smaller radius R < R^), it has no effect 
(continue to save the trajectory, but do not modify the properties 



of the object). This corresponds to a decline in the "background" 
density field - however, because the object is self-gravitating, this 
cannot simply "random walk" the collapsed region back into be- 
ing uncollapsed. By definition, gravity will prevent such expan- 
sion modulo some stronger forces applied in the model (discussed 
below). However, if the trajectory now crosses the barrier Sc{R) 
at some larger Rc.mw > R^, this corresponds to a mass growth 
event for the collapsed object. The mass of the cloud should be 
updated to Mc,new — ^ (47r/3)pc(/?c,new)^c.ncw > Mc, and the first- 
crossing/collapse radius updated to Rc — >■ Renew Unlike the case 
with dark matter halos (where all mass is locked into halos, so 
every growth event is a halo-halo merger), the fact that there is 
un-collapsed mass means that some of these events correspond to 
cloud-cloud mergers, while others coiTespond to previously "dif- 
fuse" gas reaching a self-gravitating threshold. If this distinction 
is needed, the method in |Somerville & Kolatt| (|1999|l can easily be 
generalized to decompose the mass growth AM — Mc.new — A^c into 
a "progenitor cloud" and "diffuse" mass functionp^ 

(4) Apply whatever cloud-specific physics are desired, in the 
timestep At, for the population of identified bound objects. This 
is where the essence of any semi-analytic model enters. One could 
assume clouds continue to collapse under gravity, that they form 
stars (either instantaneously, or with some efficiency in time, or 
with some association with clump-clump mergers), that they form 
molecules (based on e.g. their local column densities and SFRs), 
that they disrupt on some timescale or as some function of star 
formation/feedback properties, that they accrete "diffuse" material 
(e.g. Bondi-Hoyle accretion, which as a non-local effect is not in- 
cluded in the "growth events" in step (3)). There are obviously a 
huge range of model physics than can be included. 

One particularly interesting application of this model to bound 
clouds is to consider recursively applying the same analysis within 
each cloud, to determine the bound sub-units into which it will frag- 
ment. For a given bound cloud, if the model defines some average 
density and turbulent power spectrum (for example, assuming they 
maintain their properties at collapse, virialize and contract by dis- 
sipation, etc), then the procedure to determine the mass function 
and other properties of these "sub-clumps" is exactly identical to 
the procedure for the "parent" clumps, but with the revised or re- 
normalized density/mass/velocity properties of the "parent" clump. 

(5) Repeat steps (2)-(4), to continue to evolve the system in 
time as desired. 

We also note that despite our stated assumption of steady-state 
turbulence, it is perfectly possible to make the global parameters 
of the model (e.g. densities, masses, assumed structural properties, 
turbulent power spectrum) arbitrary functions of time and/or conse- 
quences of the explicit "cloud physics" put into the model. For ex- 
ample, allowing for continuous accretion and/or gas exhaustion to 
systematically change the normalization of the density with time, 
or allowing turbulent velocities to damp in the absence of some 
feedback from clouds to "pump" them. Likewise, it is possible to 
repeat or rescale these experiments in different "intervals" corre- 
sponding to the average properties at different radii in a galaxy disk 
(corresponding to e.g. an exponential profile) so that together, the 
Monte-Carlo ensemble can represent the properties of the entire 
galaxy disk. The only implicit assumption in the above is that these 



" To first approximation, this has the same behavior as the halo case: 
namely that the "progenitor" mass function has a similar shape to the "col- 
lapsed object" mass function, here with a similar "diffuse" mass fraction. 
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Figure 10. Fraction of the total ISM gas mass and volume in bound 
clouds, as a function of the cloud lifetime (in units of the cloud cross- 
ing time). We follow a full population of clouds thi'ough a time-dependent 
"merger/fragmentation tree" constructed as described in § |7.2| When a 
bound region collapses, we allow it to remain collapsed for a time fiifeiimc = 
i'fdyn = L/ {Rc/v,[Rc]) (?(]yn is the crossing time at the moment of becom- 
ing bound); when this time expires the mass is returned to the "diffuse" 
(non-bound) ISM. A lifetime ~ 1 — 5 (dyn gives a fraction in bound units 
consistent with the observed ISM; larger values lock all mass into bound 
units (and will over-predict the GMC MF), smaller values the opposite. 



properties evolve slowly, relative to the local mixing/equilibration 
time for the turbulence (a crossing time). 

7.3 Example: The Rate of Collapse into Bound Units and 
Constraints on Cloud Lifetimes 

It is not our intention here to present a full semi-analytic model for 
clouds. However, we briefly illustrate how such a model might be 
used with a highly simplified implementation. 

We follow steps (0)-(5) above, with the standard (dimension- 
less) parameters and p — 2 power spectrum adopted throughout 
this paper. The specification of the power spectrum and assump- 
tion of marginal stability completely specify the model, except for 
the physics applied to bound objects, step (4). 

For these bound objects, we apply a toy "zero physics" model, 
with the only goal being to see the effects of different "cloud life- 
times" on the distribution of the ISM. When an object has col- 
lapsed, we allow it to remain collapsed for a "lifetime" L/cio.ss 
where fcross is the crossing time of the cloud at the time of col- 
lapse i?c/v,(i?c) ~ 1 / \/Gpc- When this time has elapsed, we "de- 
stroy" the cloud and recycle its material, in practice by "resetting" 
the associated path (setting all AJj = for the designated path). 
The trajectory then re-grows with time according to Equation |51| 
essentially randomizing the density in a crossing time. 

For any choice of L, the mass function and mass fraction in 
bound objects will eventually converge to a steady state value (in 
practice, this requires only a couple of disk crossing times). Fig- 
ure [To] shows the resulting steady-state mass fraction in collapsed 
and bound objects, as a function of L from L <C 1 to L 3> 1. 
When L ^ 1, the mass fraction in collapsed objects is negligi- 
ble, and declines oc L at lower L (as expected for systems with a 
constant formation rate). When L ^ 1, the mass in collapsed ob- 
jects quickly saturates near unity (with an exponentially suppressed 
residual non-collapsed mass). In this regime, because clouds live 
much longer then the typical cloud-cloud merger time, the mass 
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Figure 11. Mass fraction which collapses into bound sub-units per free-fall 
time (defined as % = 0.54/y'Gpo ~ /j/(t(/!) = f2~', if Q = 1), for sy.stems 
which are marginally stable on large scales, as a function of Mach number 
(normalized at lai'ge scales, Mi, = Jv[{li)). Solid line is the analytic pre- 
diction from running a suite of models as in Figure [To] in which clouds 
remain bound (are removed from the "diffuse" medium). Given the same 
global stability, systems with higher Jv[ have larger dispersions and more 
rapidly cross the collapse barrier, but the scaling is weak. Without some 
additional physics to disrupt bound objects, mass collapses at ~ Mgas/^jy,, 
independent of the maintenance of the large-scale cascade. We compare 
with the results of high-resolution simulations of turbulence. First, ideal- 
ized forced turbulent box simulations with collapse into sink particles from 
IVazquez-Semadeni et ar])2003} (hydro) and |Padoan & Nordlund|p011) 
(hydro and MHD). For each, we select the runs which most closely match 
the "marginally stable" assumption (ctvir ~ 1 for the box). Second, full disk 
galaxy simulations with no stellar feedback from Hopkins et al. ( 20Tl} with 
"collapse rates" defined as the modeled star formation rate (which occurs in 
bound gas at densities n > lOOOcm"^, on ^ 1 pc scales). As shown therein, 
these all maintain Q ~ 1 via gravitational instabilities; A^/, is taken as the 
mass-weighted average for the disk gas averaged over a spatial scale 
= h. Error bars show the scatter in both quantities. The excursion-set model 
successfully predicts the results of fully nonlinear hydrodynamic simula- 
tions, within the scatter between different simulations/realizations. 



function also shifts to higher and higher masses (roughly shifting 
the break/maximum GMC mass Mbieak oc L). 

Only choices with L ~ 1 — 5 yield reasonable total collapsed 
masses in steady state (of order tens of percent, but with order 
tens of percent of the mass also in the inter-clump medium), and 
agreement with the observed GMC mass function. This is easy to 
understand: because the density distribution evolves on a cross- 
ing time, the rate of addition of mass to the GMC population is 
~ exp (-i'(/iraax) V2) Mgas / ?dyn (/jmax) wherc /jniax ~ h represents the 
most unstable wavelength, where v(hm:ii) is order unity. So the life- 
time for an appropriate steady-state should be L ~ exp (i'n,ax/2) ~ a 
few. 

This relates directly to idealized hydrodynamic simulations of 
turbulent boxes with self-gravity. These experiments have found 
that even when a forcing term is included to maintain the turbulent 
cascade at all times (for a box which is globally stable against col- 
lapse), a large fraction (tens of percent) of the mass in the box will 
reach densities where it becomes self-gravitating (presumably turn- 
ing into stars, if there is no feedback) in a free-fall time (see e.g. the 
discussion in |Padoan & Nordlund|201 1) 1. Here, we have calculated 
the exact same quantity analytically (on a galaxy-wide scale). 

We can estimate the rate of collapse, in the absence of feed- 
back, by simply assuming clouds are arbitrarily long-lived and then 
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calculating the time for some fraction of the mass to be bound into 
clouds. If we perform this exercise as a function of the dimension- 
less Mach number Mi, (for & p = 2 spectrum), we obtain 

^""fr ^1.5- 0.34Vln(l+3M„/8) (62) 

Where we define fconsumption as the time to 1/2 of gas consumed 
and tiyrtiji) = h/{vf(h))''^^ ~ fi"'. Note the weak and positive de- 
pendence of the collapse rate on Mh - this comes from our as- 
sumption of marginal stability for the disk as a whole - at a fixed 
stability level, larger A4 means a broader density PDF, and so in- 
creases the collapse rate. We compare the resulting collapse rate 
as a function of A4 to full numerical hydrodynamic simulations: 
both simulations of small-scale, idealized turbulent boxes (in which 
self-gravitating regions at the resolution limit are identified as sink 
particles), and large-scale simulations of galaxy disks (without stel- 
lar feedback) in which self-gravitating regions become "star parti- 
cles." In all cases we compare models with marginal stability on the 
largest scales. Our analytic calculation is in good agreement with 
the full numerical calculation. 



8 DISCUSSION 

We have used the fact that the ISM is super-sonically turbulent on a 
wide range of scales to develop a rigorous excursion-set model for 
the formation, structure, and time evolution of gas structures (e.g. 
GMCs, massive clumps/cores, and voids) in the ISM. We derive the 
conditions for self-gravitating collapse in a turbulent medium ap- 
plicable on both small scales (the Jeans condition) and large scales 
(the Toomre criterion); together with the assumption that the den- 
sity distribution in super-sonic turbulence is approximately lognor- 
mal, we use this to derive the statistical properties of the smoothed 
density field on all scales as a function of smoothing scale R. We 
show, then, that (with some appropriate modifications from the 
standard cosmological case) this becomes a well-defined barrier 
crossing problem (albeit one with a complicated barrier structure), 
for which the full methodology of excursion set theory can be ap- 
plied. 

We use this model to calculate the mass function of bound gas 
structures (over the entire dynamic range from near the sonic length 
to masses well above the Jeans mass). We do so in a rigorous man- 
ner that explicitly resolves the "cloud in cloud" problem. We show 
that this agrees extremely well with observed GMC mass functions 
in the MW and other nearby galaxies. This prediction is nearly in- 
dependent of any free parameters, with the only input being the 
mass and size of the galaxies. Even galaxies such as M33, which 
has been extensively discussed as apparently exhibiting a deviant 
GMC mass function slope, are accurately predicted. The generic 
properties of the mass function are rigorously derived: an exponen- 
tial cutoff above the Jeans mass (because large-scale density fluc- 
tuations are suppressed by disk shear) and a faint-end power-law 
slope close to, but slightly shallower than, —2 (which deviates log- 
arithmically with mass). It is near —2 (equal mass on all scales) 
generically because the collapse threshold (being relative to Inp) 
is only a logarithmic function of scale and gravity is scale free; but 
slightly shallower because collapse is more difficult on small scales 
for any realistic turbulent power spectrum. We show this is robust 
to large variations in mach number and velocity power spectrum 
shape, and even to large deviations from exact log-normality in the 
density PDF. 

The same model also predicts the linewidth-size and size-mass 



relations of these clouds, in good agreement with observations. The 
linewidth-size relation slope is a generic result of the assumed tur- 
bulent power spectrum, but its normalization is predicted by the 
assumption that the disk must be globally stable; the size-mass re- 
lation follows from the collapse criteria. Second-order corrections 
(from e.g. disk shear) make both less sensitive to the turbulent in- 
dex in the range p ~ 5/3 — 2. Residuals from these relations nat- 
urally emerge as a function of the galaxy surface density, in good 
agreement with recent comparisons of GMC properties in the MW 
outskirts and center and in high-redshift galaxies. 

The excursion set theory also allows us to rigorously predict 
the spatial correlation function and clustering properties of these 
clouds; we predict that most of the mass in clouds should be weakly 
biased (i.e. trace the overall gas density), but the most massive 
clouds will preferentially be biased towards large-scale overdensi- 
ties (e.g. spiral arms). We construct the auto-correlation function of 
GMCs from the catalogue of clouds in M33 and show this agrees 
extremely well with that predicted for clouds of the same mass. 
If we assume that young star clusters should more or less trace 
the positions of their "parent" clouds, then we can also compare 
their clustering. We show that both the star cluster-gas mass cross- 
correlation function, and the star cluster auto-correlation function 
observed in the youngest clusters in the Antennae and M5 1 agrees 
well with that predicted, over the observed range of scales R <0.lh 
toR> lOh. 

Using similar logic as applied to the GMC mass function, we 
can predict the size and mass distributions of under-dense "bub- 
bles" in the ISM. We show that a large fraction of the ISM should 
be in highly under-dense "bubbles" with p < 0.01 — 0.1 po, and 
that the characteristic size should scale with scale height h, as a 
natural consequence of turbulent fluctuations. These require no ad- 
ditional "power source" other than whatever maintains the turbu- 
lence. If we consider the distribution of bubble/hole sizes below 
a density threshold such that they should be easily ionized by the 
galactic background, we show that this agrees very well with the 
HI "hole" size distribution observed in nearby galaxies such as the 
SMC, M3 1 , and Holmberg II. The energetic cost of "creating" these 
bubbles is negligible, as compared to the nominally large PdV work 
required if they were excavated by e.g. SNe explosions, and they do 
not require any internal stars/star clusters to power their expansion. 
Even if some are powered in this manner, it is clear that many are 
not. This resolves a long-standing problem, as follow-up observa- 
tions of these "holes" have consistently failed to observe SNe rem- 
nants or other evidence of young stellar populations "powering" the 
hole expansion (see e.g. |Rhode et al.|I999||Weisz et al.|2009| l. We 
stress that turbulence alone will not explain the gas in bubbles be- 
ing hot: the fraction of holes predicted to have a cooling time much 
longer than a dynamical time from turbulent shocks alone is small. 
But it will explain their sizes, expansion, and densities. Where they 
are heated, it requires a much smaller amount of energy, at that 
point, to simply "leak into" the bubble. 

We generalize the excursion-set model of the ISM to allow 
the construction of time-dependent "merger/fragmentation trees" 
which can be used to follow the evolution of clouds and construct 
semi-analytic models for GMC and star-forming populations. We 
provide explicit recipes to construct these trees. We use a simple 
example to show that, if clouds were not destroyed by some feed- 
back process in a timescale ~ 1 — 5 crossing times, then all the ISM 
mass would be "consumed" (collapsing to arbitrarily high densi- 
ties in bound objects) even if the large scale turbulent cascade were 
maintained. Absent such feedback, we show that our analytic calcu- 
lation can predict with reasonable accuracy the collapse rates seen 
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in full nonlinear hydrodynamic (and MHD) simulations of both tur- 
bulent boxes and galaxies over a wide range of characteristic Mach 
numbers. 

It is striking that we can predict so many properties of a 
highly complex, chaotic, and - unlike the cosmological case ~ 
fully nonlinear system with a single model. This suggests that a 
wealth of properties of the ISM and GMC populations are generic 
consequences of collapse in a supersonically turbulent medium 
with a characteristic "scale" set by gravitational instability in a 
gaseous disk. This explains why different simulations fOstrikerl 
erar]|200Tl |Dobbs||2508l |Boumaud et~aL] [2007; Tasker & Tan 



2009[ Tasker| 201 1 1 have been able to reproduce various aspects 
of these observations, despite including very different physics for 
cooling/feedback/star formation/magnetic fields, and in some cases 
clearly failing to reproduce other (probably feedback-dependent) 
properties such as the observed star formation rates and galactic 
winds (see e.g. the discussion in [Hopkins et al.||2012[ and refer- 
ences therein). What is remarkable is that our theory allows us to 
calculate these nonlinear properties analytically, over a large dy- 
namic range, and in quantitative agreement with the observations. 

We should also stress that this model does not necessarily im- 
ply or require that the ISM structure be rigorously self-similar or 
fractal: that may be true, but it is a much stronger statement about 
the structure of S{R) and B[R) (compared to our assumptions). In 
our default model those happen to be approximately scale free over 
some range, but there are at least two scales - the sonic length and 
disk scale height - above/below which behavior changes. The ap- 
plication of this model also does not necessarily imply the ISM 
is "hierarchical" either in the cosmological sense or the sense of 
IVazquez-Semademl ( 1994 1 (see § [TJ. In fact in the cosmological 
sense of the term, the predicted structure is more appropriately 
"anti-hierarchical," in that collapse tends to proceed "top-down." 
Large scales ~ h are most unstable and contain most of the turbu- 
lent power, and we have shown that most self-gravitating objects 
on small scales are formed by "fragmenting out" of larger struc- 
tures (i.e. form within parent GMCs; these are the low-mass GMCs 
predicted if we ignore the cloud-in-cloud problem, which are much 
more abundant than isolated counterparts). In contrast, in the cos- 
mological case, small structures form first, and "subhalos" are only 
a small fraction of the population at most masses. 

There are many interesting potential future directions for these 
models. Many of our assumptions can be made more general, and 
the model made more accurate. For example, if the gas is not 
isothermal, or when magnetic fields and global gravitational forces 
are strong, some deviations from lognormality are expected. We 
have argued that should not qualitatively change our conclusions, 
but it is possible to rigorously treat this regime, by extending the ex- 
cursion set formalism to non-Gaussian fields (as developed in e.g. 



Matarrese et al.|2000[[Afshordi & Tolley|2008[|Maggiore & Riotto| 



2010b I. The Monte Carlo excursion set approach is also completely 
generalizeable to treat correlated random walks, so that arbitrary 
higher-order structure functions (i.e. correlations between fluctua- 
tions on different scales) can be incorporated - it is only a con- 
venient simplifying assumption to assume strict locality (see the 
review in Zentner 2007 1. Near and below the sonic length (or in 
the warm/diffuse ISM), when the turbulence is sub/transsonic, ad- 
ditional corrections to the power spectrum could be included. Mag- 
netic fields can also be included as more than just a correction to 
the effective sound speed, if their power spectrum is well deter- 
mined (see e.g. |Kim et al.|2002^ . We have also assumed that the 
density and velocity field are not directly coupled, but it is in prin- 
ciple straightforward to allow for a direct correlation between the 



local density and velocity fluctuations, to follow both with a higher- 
dimensional random walk, and to incorporate this in the collapse 
criterion (see e.g. |Sheth & Tormen|2002 l. We have neglected large- 
scale gradients in e.g. the disk surface density profile; this should 
not be important for most GMCs since their sizes are < h, but a 
more rigorous calculation of global properties (e.g. the large-scale 
spatial distribution of clouds) could consider each radial annulus 
in turn with some global mass profile. Our derivation of the col- 
lapse barrier also assumes spherical collapse, when in fact most 
GMCs are ellipsoidal or triaxial. In the cosmological context, ellip- 
soidal collapse is a well-studied problem (Sheth et al. 2001 1, and 
can be incorporated via an appropriate change of the barrier shape 
(although the appropriate parameters are usually determined by ref- 
erence to numerical simulations); however, if the cosmological case 
is any guide, the differences should be small (tens of percent level). 

The lack of dependence of many of the predicted observables 
on the detailed properties of turbulence is, in one sense, reassuring 
(and explains agreement between previous models with different 
physics). On the other hand, unfortunately, it means that observa- 
tions have a limited ability to discriminate between these differ- 
ent physical scenarios. Of course, the ISM is not all highly super- 
sonic, and there will be some regimes in which the model here is 
not appropriate. Implicitly, our model assumes that the ISM can 
cool efficiently (cooling time short relative to the dynamical time), 
so that the turbulence remains super-sonic. This is easily satisfied 
inside the radii that include most star formation in galactic disks. 
However, at sufficiently large radii and very low gas densities, the 
galactic and cosmological ionizing backgrounds maintain gas disks 
as fully ionized with Q> \. Even in the star-forming disk, there are 
bubbles of hot gas that may escape before cooling. And a signifi- 
cant fraction of the mass in the "warm" ISM medium is turbulent 
and bound, but has comparable thermal sound speeds and turbulent 
velocities (i.e. is transsonic rather than super-sonic). In this regime, 
it is necessary to account explicitly for the effects of heating and 
cooling, for example as in the model of |Ostriker et al.| |2010 1, since 
only regions pushed to a critical density where cooling becomes ef- 
ficient will behave as we assume. Even in this regime, however, the 
internal structure of those cooling regions/GMCs should be treat- 
able with the model here, but it should be emphasized that this 
model is most applicable to the cold/rapidly cooling gas rather than 
the extended low-density gas. 

There is also a huge space of physical predictions which can 
be explored with this model that we have not yet addressed, some 
of which may be more sensitive tests of the properties of turbu- 
lence and ISM structure formation. Using the time-dependent for- 
mulation, the growth of of GMCs via turbulent density fluctuations, 
mergers, and accretion can be rigorously analytically calculated. 
By allowing for global evolution of self-gravitating regions, it is 
possible to self-consistently follow features that nominally appear 
to contradict the model - for example, following a simple model 
whereby, once "detached" from the background, a cloud collapses 
spherically will naturally predict a power-law tail to high densi- 
ties (in collapsing regions), even though we can continue to use 
the same treatment for each cloud internally. Together with a semi- 
analytic model for star formation in such units, their destruction via 
feedback can also be followed analytically in a self-consistent sta- 
tistical model for the population. Global feedback effects can also 
be predicted - for example, many authors have used the cosmo- 
logical formulation of the problem to study the reionization history 
of the Universe and evolving size distribution of HII regions (e.g. 
|Haiman et al.|2000[[Furlanetto et al.|2004^ . It is straightforward to 
adapt their approach to the problem here to predict the properties 
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of galactic HII regions, SNe blastwaves, and ionizing photon es- 
cape fractions. The model can be extended iteratively (downwards 
in scale) within GMCs, to calculate the properties of dense collaps- 
ing subregions (cores). Extended sufficiently in scale, the model 
can even be used to predict the stellar IMF in each subregion, fol- 
lowing [Hennebell£|^|Ch^ri^ - with a model determined 
on galactic scales. These scales, being closer to the sonic length, 
should exhibit much stronger dependence on the actual turbulent 
structure than the galactic-scale quantities we calculate here (as 
seen in other analytic calculations and simulations; see |Ballesteros-| 
[Paredes et a"Ll[20g6l [Hennebelle & Chabrier|[2009) , and might be 
used to break degeneracies between different models for the ISM 
microphysics. 
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